IIR滤波器设计

说明 1：由于用Matlab设计出系数之后，要在实验板上面跑，所以取两个正弦信号。

       Fs=1000000；               %采样率是1M

       N=1024;                    %设置成1024点FFT

       x1=sin(2*pi*22500*t);      %有时候作为信号  22500Hz
       x2=sin(2*pi*45000*t);     %有时候作为噪声  45000Hz

    2：主要设计了巴特沃斯，切比雪夫I型，椭圆滤波器的高通，低通，带通，带阻。

       并将高通和低通设计成固定的两阶，带通和带阻设置成一阶，需要高阶的话可

       以级联。当然也可以按照模拟滤波器的函数，计算出需要的截止频率和需要的 

       阶数，但是效果不是特别的好，所以就不采用这种方法了。

    3. 滤波后的效果，主要的通过查看输出后的波形和信噪比得出滤波的效果好

      坏。

    4. 以巴特沃斯作为说明，其它的滤波器的设计类似。

    5. 这里的信噪比SNR的计算

       实际的应用中信噪比的计算式，以放大器信噪比为例

       给放大器一个标准信号，通常是0.775Vrms或2Vp-p@1kHz,调整放大器的放大倍数使其达到最大不失真

       输出功率或幅度（失真的范围由厂家决定,通常是10%，也有 1%），记下此时放大器的输出幅Vs，然后

       撤除输入信号，测量此时出现在输出端的噪声电压，记为Vn，再根据SNR=20LG(Vn/Vs)就可以计算出信

       噪比了。Ps和Pn分别是信号和噪声的有效功率，根据SNR=10LG（Ps/Pn）也可以计算出信号比。

       我采用的信噪比的计算方法是

       Ps=sum(x1.^2);                   %信号的总功率  X1是信号
     Pu=sum((y2-x1).^2);              %剩余噪声的功率  y2是滤波后的信号
     SNR=10*log10(Ps/Pu);             %可以认为是平均的功率得到的信噪比

     还有一个特别的需要注意的是计算差分方程不要使用这个函数filter，而是得使

     用filtfilt 这样的话计算出的信噪比是正常的。 

一 巴特沃斯低通滤波器

   

fs=1000000;                    %设置采样频率 1k
N=1024;
n=0:N-1;
t=0:1/fs:1;
f=n*fs/N;
x1=sin(2*pi*22500*t);          %信号
x2=sin(2*pi*45000*t);          %噪声
x=x1+x2;
y=fft(x,N);
subplot(221);
plot(t,x);title('原始信号');grid on;
subplot(222);
plot(f,abs(y));title('原始信号FFT');

Wc=2*30000/fs;                               

[b,a]=butter(2,Wc);
% y2=filter(b,a,x);

y2=filtfilt(b,a,x);
Ps=sum(x1.^2);                    %信号的总功率
Pu=sum((y2-x1).^2);               %剩余噪声的功率
SNR=10*log10(Ps/Pu);             %可以认为是平均的功率得到的信噪比

y3=fft(y2,N);
subplot(223);                     %2阶IIR滤波后的FFT
plot(f,abs(y3));xlabel('频率/Hz');
ylabel('振幅');title('滤波后信号FFT');
[H,F]=freqz(b,a,512);
subplot(224);
plot(F/pi,abs(H));xlabel('归一化频率');        %绘制绝对幅频响应
ylabel('幅度');
title(['Order=',int2str(2),'    SNR=',num2str(SNR)]);grid on;
set(gcf, 'color', 'white');

二，巴特沃斯带通滤波器的设计

fs=1000000;                         %设置采样频率 1k
N=1024;
n=0:N-1;
t=0:1/fs:1;
f=n*fs/N;
x1=sin(2*pi*22500*t);                %这里X1是噪声  需要通过带阻滤波滤除
x2=sin(2*pi*45000*t);                %这里X2是信号，
x=x1+x2;
y=fft(x,N);
subplot(221)
plot(t,x);title('原始信号');grid on;
subplot(222);
plot(f,abs(y));title('原始信号FFT');
                                     %方法一
%  Ws= [15000*2 30000*2]/fs;         %Ws>Wp的时候是带阻，Ws<Wp的时候是高通，      
%  Wp=[10000*2   35000*2]/fs;        %这里是通过计算函数buttord计算的最小的阶数和截止频率      
%  [n,Wn]=buttord(Wp,Ws,3,20);       %计算出来的效果不怎么的好
%  [b,a]=butter(n,Wn);

Wn=[43000*2 47000*2]/fs;            %方法二  直接的锁定带阻截止频率的范围
[b,a]=butter(1,Wn);

% y2=filter(b,a,x);                  %用filtfilt更加的有利于计算信噪

y2=filtfilt(b,a,x);
Ps=sum(x2.^2);                       %信号的总功率
Pu=sum((y2-x2).^2);                  %剩余噪声的功率
SNR=10*log10(Ps/Pu);                 %可以认为是平均的功率得到的信噪比

y3=fft(y2,N);
subplot(223);                        %2阶IIR滤波后的FFT
plot(f,abs(y3));
xlabel('频率/Hz');
ylabel('振幅');
title('滤波后信号FFT');grid on;
[H,F]=freqz(b,a,512);
subplot(224);
plot(F/pi,abs(H));
xlabel('归一化频率');               %绘制绝对幅频响应
ylabel('幅度');
title(['Order=',int2str(1),'    SNR=',num2str(SNR)]);grid on;
set(gcf, 'color', 'white');
% figure;
% freqz(b,a,512);
% figure;
% plot(t,y2);

三，巴特沃斯带阻滤波器

fs=1000000;                               %设置采样频率 1k
N=1024;
n=0:N-1;
t=0:1/fs:1;
f=n*fs/N;
x1=sin(2*pi*22500*t);                     %这里X1是噪声  需要通过带阻滤波滤除
x2=sin(2*pi*45000*t);                     %这里X2是信号，
x=x1+x2;
y=fft(x,N);
subplot(221)
plot(t,x);title('原始信号');
subplot(222);
plot(f,abs(y));title('原始信号FFT');
                                 %方法一
%  Ws= [15000*2 30000*2]/fs;     %Ws>Wp的时候是带阻，Ws<Wp的时候是高通，      
%  Wp=[10000*2   35000*2]/fs;   %这里是通过计算函数buttord计算的最小的阶数和截止频率      
%  [n,Wn]=buttord(Wp,Ws,3,20);  %计算出来的效果不怎么的好
%  [b,a]=butter(n,Wn);

Wn=[20000*2 25000*2]/fs;                 %方法二  直接的锁定带阻截止频率的范围
[b,a]=butter(2,Wn,'stop');

% y2=filter(b,a,x);                       %用filtfilt更加的有利于计算信噪

y2=filtfilt(b,a,x);
Ps=sum(x2.^2);                            %信号的总功率
Pu=sum((y2-x2).^2);                       %剩余噪声的功率
SNR=10*log10(Ps/Pu);                      %可以认为是平均的功率得到的信噪比

y3=fft(y2,N);
subplot(223);                             %2阶IIR滤波后的FFT
plot(f,abs(y3));
xlabel('频率/Hz');
ylabel('振幅');
title('滤波后信号FFT');grid on;
[H,F]=freqz(b,a,512);
subplot(224);
plot(F/pi,abs(H));
xlabel('归一化频率');                     %绘制绝对幅频响应
ylabel('幅度');

title(['Order=',int2str(2),'    SNR=',num2str(SNR)]);grid on;
set(gcf, 'color', 'white');

四，巴特沃斯高通滤波器

 

fs=1000000;                          %设置采样频率 1k
N=1024;
n=0:N-1;
t=0:1/fs:1;
f=n*fs/N;
x1=sin(2*pi*22500*t);               %噪声
x2=sin(2*pi*45000*t);               %信号
x=x1+x2;
y=fft(x,N);
subplot(221)
plot(t,x);title('原始信号');
subplot(222);
plot(f,abs(y));title('原始信号FFT');

%  Ws=10000*2/fs
%  Wp=30000*2/fs;                                    

%  [n,Wn]=buttord(Wp,Ws,3,20);
%  [b,a]=butter(n,Wn,'high');
[b,a]=butter(2,32000*2/fs,'high');

% y2=filter(b,a,x);

y2=filtfilt(b,a,x);
Ps=sum(x2.^2);                               %信号的总功率
Pu=sum((y2-x2).^2);                          %剩余噪声的功率
SNR=10*log10(Ps/Pu);                         %可以认为是平均的功率得到的信噪比

y3=fft(y2,N);
subplot(223);                                %2阶IIR滤波后的FFT
plot(f,abs(y3));xlabel('频率/Hz');
ylabel('振幅');title('滤波后信号FFT');grid on;
[H,F]=freqz(b,a,512);
subplot(224);
plot(F/pi,abs(H));xlabel('归一化频率');          %绘制绝对幅频响应
ylabel('幅度');title(['Order=',int2str(2),'    SNR=',num2str(SNR)]);grid on;
set(gcf, 'color', 'white');
% figure;
% % freqz(b,a,512);
% plot(t,y2);