Matlab中，编写一个用模拟低通滤波器和切比雪夫模拟低通滤波器对加噪之后的原信号处理的程序

clc;clear all;

%归一化模拟切比雪夫I型低通滤波器的设计

Wp=2*pi*1000;Ws=2*pi*1500;rp=3;rs=30;%设计滤波器的参数    

wp=1;ws=Ws/Wp;                       %频带变换得到归一化滤波器

[N,wc]=cheb1ord(wp,ws,rp,rs,'s');    %计算滤波器阶数和3dB截止频率  

[z,p,k]=cheb1ap(N,rp);               %计算归一化滤波器的零极点

[b,a]=zp2tf(z,p,k);                  %计算归一化滤波器系统函数的系数

w0=0:0.05*pi:2*pi;

[h0,w0]=freqs(b,a,w0);               %求频率响应 

figure(1)

plot(w0,20*log10(abs(h0)),'k');

title('归一化模拟且比雪夫I型低通滤波器');

xlabel('频率f/Hz');ylabel('幅度/dB');grid;

%一般低通切比雪夫低通滤波器的设计

[B,A]=lp2lp(b,a,Wp);                 %将归一化滤波器转换为一般模拟滤波器

w1=0:2*pi*100:2*pi*30000;

[h1,w1]=freqs(B,A,w1);

figure(2)

plot(w1/(2*pi),20*log10(abs(h1)),'k');

title('一般模拟且比雪夫I型低通滤波器');

xlabel('频率f/Hz');ylabel('幅度/dB');grid;

%冲激响应不变法设计低通巴特沃斯数字滤波器

Fs=10000;                            %采样频率 

Wp1=Wp/Fs;Ws1=Ws/Fs;

rp1=3;rs1=30;                        %设计滤波器的参数  

[N1,Wn]=cheb1ord(Wp,Ws,rp1,rs1,'s')  %计算滤波器阶数

[b1,a1]=cheby1(N1,rp1,Wn,'s');       %样本AF的系数函数的分子分母系数

[bz,az]=impinvar(b1,a1,Fs);          %冲激响应不变法从AF到DF变换

[C1,B1,A1]=dir2par(bz,az)            %直接形式转换为并联型

w2=[Wp1,Ws1];                        %数字临界频率

[H,f]=freqz(bz,az);                  %绘制数字滤波器的幅频特性和相频特性

[db,mag,pha,grd,f]=freqz_m(bz,az);   %扩展函数检验滤波器的其他特性

figure(3)

plot(f/pi,db,'k');

title('冲激响应不变法设计低通切比雪夫I型数字滤波器');

xlabel('角频率w/pi');ylabel('振幅/dB');

axis([0,0.35,-30,5]);grid;

%用设计好的滤波器对信号进滤波处理

figure(4)

f1=500;f2=4000;                      %输入信号的频率

N=100;                               %数据长度

dt=1/Fs;n=0:N-1;t=n*dt;              %采样间隔和时间序列

x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号

subplot(2,1,1),plot(t,x),title('输入信号');

y=filtfilt(bz,az,x);                 %用滤波器进行滤波处理

y1=filter(bz,az,x);                  %进行滤波处理

subplot(2,1,2),plot(t,y,t,y1,':'),title('输出信号');xlabel('时间/s');

legend('filtfilt','filter')          %加图例

freqz_m.m文件

function[db,mag,pha,grd,w]=freqz_m(b,a)

[H,w]=freqz(b,a,1000,'whole');

mag=abs(H);

db=20*log10((mag+eps)/max(mag));

pha=angle(H);

grd=grpdelay(pha);

dir2par.m文件

function [C,B,A] = dir2par(b,a)

M=length(b);N=length(a);

[r1,p1,C]=residuez(b,a);

p=cplxpair(p1,10000000*eps);

I=cplxcomp(p1,p);

r=r1(I);

K=floor(N/2);B=zeros(K,2);A=zeros(K,3);

if K*2==N;

    for i=1:2:N-2

        Brow=r(i:1:i+1,:);

        Arow=p(i:1:i+1,:);

        [Brow,Arow]=residuez(Brow,Arow,[]);

        B(fix((i+1)/2),:)=real(Brow);

        A(fix((i+1)/2),:)=real(Arow);

    end

    [Brow,Arow]=residuez(r(N-1),p(N-1),[]);

    B(K,:)=[real(Brow) 0];A(K,:)=[real(Arow) 0];

else

    for i=1:2:N-1

        Brow=r(i:1:i+1,:);

        Arow=p(i:1:i+1,:);

        [Brow,Arow]=residuez(Brow,Arow,[]);

        B(fix((i+1)/2),:)=real(Brow);

        A(fix((i+1)/2),:)=real(Arow);

    end

end

cplxcomp.m文件

function I=cplxcomp(p1,p2)

I=[];

for j=1:1:length(p2)

    for i=1:1:length(p1)

        if(abs(p1(i)-p2(j))<0.0001)

            I=[I,i];

        end

    end

end

I=I';