自适应IIR格型滤波器的Matlab实现

%IIR算法的Matlab实现，作者Roy,时间：2022-05-11 %本文将以正弦信号为例，在正弦信号上加高斯噪声和其他频段的正弦信号干扰 %通过rife算法估计出干扰频率，采用自适应滤波器IIR算法，滤除正弦波干扰信号。 clc,clear,close all; f=200; %信号频率为 A=4; %设置正弦波幅值 N=1024; %采样点个数 n=0:N-1; %采样序列 fs=2000; %采样频率为 fs=2K ts=1/fs; %采样间隔 T=N*ts; %采样时间 m=-fs/2:(fs/N):(fs/2-(fs/N)); %其中fs/N为采样分辨率 % % % f1=20; s_n=sin(2*pi*f1*n*ts); subplot(4,2,1);plot(n*ts,s_n); title('原始信号：sin(2*pi*f1*n*ts)'); subplot(4,2,2);plot(m,abs(fftshift(fft(s_n)))); title('频谱分析'); SNR=10; s_n=awgn(s_n,SNR); %SNR参数，可参考awgn函数 subplot(4,2,3);plot(n*ts,s_n); title('原始信号+高斯噪声'); subplot(4,2,4);plot(m,abs(fftshift(fft(s_n)))); title('频谱分析'); i_n=A*sin(2*pi*f*n*ts); %生成正弦波信号,干扰信号 xn=s_n+i_n; %添加干扰信号 subplot(4,2,5);plot(n*ts,xn); title('原始信号+高斯噪声+窄带干扰') subplot(4,2,6);plot(m,abs(fftshift(fft(xn)))); title('频谱分析'); %%%使用rife算法进行频率估计 yw=abs(fft(xn)); %傅里叶变换 % subplot(4,1,4);plot(n,yw); [valueMax,posMax] = max(yw(1:N)); %得到谱线最大值与坐标 if( yw(posMax+1) >= yw(posMax-1)) %判断r的值 r=1; else r=-1; end delta = yw(posMax+r)/(valueMax+yw(posMax+r)); fc = 1/T*(posMax-1+r*delta); %输出估计的频率值fc,一开始没有减1，导致误差比较大，可以将1/T*(posMax-1)、1/T*(posMax)分别输出与图像对比 result = strcat('Rife算法估计出的信号频率为',num2str(fc),'Hz；','实际信号频率是：',num2str(f),'Hz；'); disp(result); % % % % % % % % % IIR滤波器的算法实现 beta=cos(2*pi*fc/fs); alpha = 0.1; mul_1=0; mul_2=0; mul_3=0; add_1=0; add_2=0; add_3=0; add_4=0; add_5=0; add_6=0; lag_1=0; lag_2=0; for nn = 1 : N mul_1=xn(nn)/2; add_1=mul_1-lag_2; mul_2=alpha*add_1; add_2=mul_2+mul_1; add_3=add_2-lag_1; mul_3=(-beta)*add_3; add_4=add_2+mul_3; add_6=mul_2+lag_2; y1(nn)=add_6+mul_1; %%%%其中数据的滞后通过将变量滞后一轮输出实现 %%%%保证数据的滞后使用的数据是上一轮的，这三行代码顺序很重要 add_5=lag_1+mul_3; lag_2=add_5; lag_1=add_4; end % b1 =0; % b0 = 0; % f1 = 0; % f0 = 0; % adder1 = 0; % multer1 = 0; % adder2 = 0; % multer2 = 0; % adderb0 = 0; % y1=iir(xn,fc,fs,100); 测试封装函数 % for mm = 1 : N % xtemp = xn( mm )/2; % b1 = b0 + multer2; % b0 = f0; % adder1 = xtemp - b1; % multer1 = alpha * adder1; % f1 = xtemp + multer1; % adder2 = f1 - b0; % multer2 = (-beta) * adder2; % f0 = f1 + multer2; % y1(mm) = b1 + multer1 + xtemp; % end subplot(4,2,7);plot(n*ts,y1); title('滤除窄带干扰后的信号') subplot(4,2,8);plot(m,abs(fftshift(fft(y1)))); title('频谱分析');