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*230000*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*247000*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*230000*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*225000*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);