ddt.m

几种数字滤波器

%1． 利用傅立叶级数展开的方法，自由生成所需的x(t)；
%2．	通过选择不同的采样间隔T（分别选T>或<1/2fc），从x(t)获得相应的x(n)（作出x(n)图形）；
%3．	对获得的不同x(n)分别作傅立叶变换，分析其频率响应特性（给出幅频与相频特性曲线）；
%4．	利用巴特沃思、切比雪夫或椭圆滤波器设计数字滤波器（滤波特性自定），要求通过改变滤波器参数或特性（低通、高通、带通或带阻）设计至少两种数字滤波器，分析所设计滤波器（画出频率特性曲线），并对上述给出的不同x(n)分别进行滤波（画出滤波结果），然后加以讨论；
%5．	利用窗函数设计法或频率采样法设计数字滤波器（滤波特性自定），要求通过改变滤波器参数或特性（低通、高通、带通或带阻等）设计至少两种数字滤波器，分析所设计滤波器（画出频率特性曲线），并对上述给出的不同x(n)分别进行滤波（画出滤波结果），然后加以讨论。
% 1.生成所需的x(t)；
f = 40;
t = linspace(0,50,2000);
w = 2*pi*f;
y1=2*sin(w*t);
y2=3*sin(3*w*t); %正弦叠加分量
figure(1)
subplot(3,1,1); 
plot(t,y1) ,grid on;
subplot(3,1,2);
plot(t,y2) ,grid on;
x=y1+y2;        %正弦相叠加
subplot(3,1,3);
plot(t,x) ,grid on;
xlabel('t');
ylabel('x(t)');
title('两分量叠加后的原始模拟信号');
%%% 2.通过选择不同的采样间隔T（分别选T>或<1/2fc），从x(t)获得相应的x(n)（作出x(n)图形）；
%2.1对x(t)采样 采样频率大于2fs=240 即T<1/2*fc 分析其频率响应特性（给出幅频与相频特性曲线）
%(1)采样
figure(2);
n=0:119;
T = 1/250;
xn=2*sin(w*n*T)+3*sin(3*w*n*T)
subplot(3,1,1);
stem(n,xn,'.');
xlabel('n');
ylabel('x(n)');
title('采样后的信号图，采样间隔T<1/2*fc');
%(2)快速弗利叶变换
yk = fft(xn,120);
i=0:119;
subplot(3,1,2);
stem(i,abs(yk),'.');
xlabel('(N =120 wk = 2*pi*k/N) k');
ylabel('|Xn(K)|');
title('FFT1');
%(3)相频特性
ang = angle(yk);
subplot(3,1,3);
stem(i,ang,'.');
xlabel('k');
ylabel('P(w)');
title('相频特性曲线1');
%2.2对x(t)采样 采样频率大于2fs=240 即T>1/2*fc 分析其频率响应特性（给出幅频与相频特性曲线）
%(1)采样
figure(3);
subplot(3,1,1);
T1=1/180
xn1=2*sin(w*n*T1)+3*sin(3*w*n*T1)
stem(n,xn1,'.');
xlabel('n');
ylabel('x(n1)');
title('采样后的信号图，采样间隔T>1/2*fc');
%(2)快速弗利叶变换
yb = fft(xn1,120);
i=0:119;

subplot(3,1,2);
stem(i,abs(yb),'.');
xlabel('(N = 120 wk = 2*pi*k/N) k');
ylabel('|Xb(K)|');
title('FFT2');
%(3)相频特性
ang1 = angle(yb);
subplot(3,1,3);
stem(i,ang1,'.');
xlabel('k');
ylabel('P(w)');
title('相频特性曲线2');

%4．利用巴特沃思、切比雪夫或椭圆滤波器设计数字滤波器（滤波特性自定），
%要求通过改变滤波器参数或特性（低通、高通、带通或带阻）设计至少两种数字滤波器，
%分析所设计滤波器（画出频率特性曲线），并对上述给出的不同x(n)分别进行滤波（画出滤波结果），然后加以讨论；
%%巴特沃思低通滤波器的设计
figure(4);
f=40
fs=250
fc=fs/2
Wp = f/fc;
Ws = 2*f/fc;
[n1,Wn] = buttord(Wp,Ws,3,30);
[b,a] = butter(n1,Wn);
[h,w]=freqz(b,a,512,fs); 
title('巴特沃思低通滤波器辐频特性');
subplot(3,1,1);
plot(w,abs(h));

y= filter(b,a,xn)
title('对xn进行滤波');
subplot(3,1,2);
plot(n,y) ,grid;


y1 = fft(y,120);
i=0:119;
subplot(3,1,3);
stem(i,abs(y1),'.');
xlabel('(N = 120 wk = 2*pi*k/N) k');
ylabel('|Xa(K)|');
title('滤波结果');

%巴特沃思带通滤波器的设计
figure(5 );
Wp =[2*f/fs 3*f/fs];
Ws =[1.5*f/fs,5*f/fs];
[n1,Wn] = buttord(Wp,Ws,3,120);
[b,a] = butter(n1,Wn);
[h2,w2]=freqz(b,a,512,fs); 
title('巴特沃思带通滤波器辐频特性');
subplot(3,1,1);
plot(w2,abs(h2));

y = filter(b,a,xn)
title('对xn进行滤波')
subplot(3,1,2);
plot(n,y);

y1 = fft(y,120);
i=0:119;
subplot(3,1,3);
stem(i,abs(y1),'.');
xlabel('(N = 120 wk = 2*pi*k/N) k');
ylabel('|Xa(K)|');
title('滤波结果');

%5．	利用窗函数设计法或频率采样法设计数字滤波器（滤波特性自定），
%要求通过改变滤波器参数或特性（低通、高通、带通或带阻等）设计至少两种数字滤波器，
%分析所设计滤波器（画出频率特性曲线），并对上述给出的不同x(n)分别进行滤波（画出滤波结果），然后加以讨论。
%%用海明窗设计高通滤波器 
figure(6 );
m = 20
wp1  =5.5*f/fs;
[b1 a1]=fir1(m,wp1,'high',hamming(m+1));
[h1 w1]=freqz(b1,a1,512,fs);
subplot(3,1,1);
plot(w1,abs(h1));
title('海明窗高通数字滤波器');

y = filter(b1,a1,xn)
subplot(3,1,2);
plot(n,y) ,grid on;
title('高通滤波后的信号');

y1 = fft(y,120);
i=0:119;
subplot(3,1,3);
stem(i,abs(y1),'.');
xlabel('(N = 120 wk = 2*pi*k/N) k');
ylabel('|Xa(K)|');
title('高通滤波后的信号频谱');

%%用矩形窗设计带通滤波器
figure(7);
m = 20;
wp2  =[2*f  4*f]/fs;
[b2 a2]=fir1(m,wp2,boxcar(m+1));
[h2 w2]=freqz(b2,a2,512,fs);
subplot(3,1,1);
plot(w2,abs(h2));
title('矩形窗带通数字滤波器');

y = filter(b2,a2,xn)
subplot(3,1,2);
plot(n,y) ,grid on;
title('带通滤波后的信号');

y2 = fft(y);
i=0:119;
subplot(3,1,3);
stem(i,abs(y2),'.');
xlabel('(N = 120 wk = 2*pi*k/N) k');
ylabel('|Xa(K)|');
title('带通滤波后的信号频谱');