chenyatu2004042080.m

数字信号处理的一个课程设计

wp=input('Input 通带截止频率 wp=');
ws=input('Input 阻带截止频率 ws=');
d=input('Input 通带最大波动 d=');
k=input('Input 阻带衰减 k= ');



if wp>=ws|wp>=1|ws>=1
    errordlg('wp或ws不符合要求，请重新输入，并使wp<ws&wp<1&ws<1');
    wp=input('Input 通带截止频率 wp=');
    ws=input('Input 阻带截止频率 ws=');
end  
wp=wp*pi;
ws=ws*pi;
if k<13|k>100
errordlg('k 不符合要求，请再输入，使k>=13&k>=100');
k=input('Input 阻带衰减 k=');
end
if d>3|d<0
errordlg('d 不符合要求，input d again,d<3&d>0');
d=input('Input 最大波动 d=');
end
p=[1 2 2 4];
w=ws-wp;
wc=(ws+wp)/2+eps;
if k>=13&k<31
N=ceil(p(1)*4*pi/w+eps);
elseif k>=31&k<41
N=ceil(p(2)*4*pi/w+eps);
elseif k>=41&k<57
N=ceil(p(3)*4*pi/w+eps);
elseif k>=57
N=ceil(p(4)*4*pi/w+eps);
end
if rem(N,2)==0
    N=N+1;
end
M=(N+1)/2;
w1=boxcar(N);
w2=hanning(N);
w3=hamming(N);
w4=blackman(N);

n=1:1:N;
hd=sin(wc*(n-M))./((pi*(n-M)+eps));
hd(M)=wc/pi;
h1=hd.*rot90(w1);
h2=hd.*rot90(w2);
h3=hd.*rot90(w3);
h4=hd.*rot90(w4);

[mag1,rad]=freqz(h1,1,512);
[mag2,rad]=freqz(h2,1,512);
[mag3,rad]=freqz(h3,1,512)
[mag4,rad]=freqz(h4,1,512);

a1=angle(mag1);
a2=angle(mag2);
a3=angle(mag3);
a4=angle(mag4);

figure(1)
plot(rad,20*log10(abs(mag1)))
title('矩形窗')
xlabel('w')
ylabel('对数幅值(db)')
axis tight
grid on
figure(2)
plot(rad,20*log10(abs(mag2)))
axis tight
title('汉宁窗')
xlabel('w')
ylabel('对数幅值(db)')
grid on
figure(3)
plot(rad,20*log10(abs(mag3)))
axis tight
title('海明窗')
xlabel('w')
ylabel('对数幅值(db)')
grid on
figure(4)
plot(rad,20*log10(abs(mag4)))
axis tight
title('布莱克曼窗')
xlabel('w')
ylabel('对数幅值(db)')
grid on
if  k>57
BETA=input('Input kaiser BETA=');
w5=kaiser(N,BETA);
h5=hd.*rot90(w5);
[mag5,rad]=freqz(h5,1,512);
a5=angle(mag5);
figure(5)
plot(rad,20*log10(abs(mag5)))
title('凯塞窗')
xlabel('w')
ylabel('对数幅值(db)')
vaxis([0 pi -160 5]);
grid on
end