fig731-2.m

高西全2009年最新版<数字信号处理>源码程序

%《数字信号处理--原理、实现及应用》第七章 图7.3.1 程序：fig731.m
% 电子工业出版社出版  高西全  丁玉美  阔永红合著 2006年8月
% 用频率采样法设计FIR低通滤波器
clear;close all
%============================================================================================
%采样点数N=13
N=15;wc=pi/3;
Np=fix(wc/(2*pi/N));Ns=N-2*Np-1; %Np+1为通带[0,wc]上采样点数,Ns为阻带[wc,2*pi-wc]上采样点数
Ak=[ones(1,Np+1),zeros(1,Ns),ones(1,Np)];     %幅度采样向量A(k)
thetak=-pi*(N-1)*(0:N-1)/N;                 %相位采样向量θ(k) 
Hk=Ak.*exp(j*thetak);    %构造频域采样向量H(k)
hn=real(ifft(Hk));      %h(n)=IDFT[H(k)]
%==============================================================================================
%绘图7.3.1
Aw=[1,1,0,0,1,1];wk1=[0,wc,wc,2*pi-wc,2*pi-wc,2*pi]/pi;
subplot(3,2,1);plot(wk1,Aw);axis([0,2,-0.2,1.2])
xlabel('ω/π');ylabel('H_g(k)');title('(a)')
hold on
wk2=[0:N-1]*2/N;
plot(wk2,Ak,'.');
n=0:N-1;
subplot(3,2,2);stem(n,hn,'.');axis([0,15,-0.2,0.5])
xlabel('n');ylabel('h(n)');title('(b)')
Hw=fft(hn,1024);
wk=2*pi*[0:1023]/1024;
Hgw=Hw.*exp(j*wk*(N-1)/2);
%=====================================================
%计算N=15时的通带最大衰减rp15和阻带最小衰减rs15
rp15=max(20*log10(abs(Hgw)))
hgmin=min(real(Hgw))
rs15=20*log10(abs(hgmin))
%========================================================
subplot(3,2,5);plot(wk/pi,real(Hgw));axis([0,2,-0.2,1.2])
hold on;
plot(wk1,Aw,':');plot(wk2,Ak,'.');
xlabel('ω/π');ylabel('Hg(ω)');title('(c)')
subplot(3,2,6);plot(wk/pi,20*log10(abs(Hgw)));axis([0,1,-30,3]);grid on
xlabel('ω/π');ylabel('20lg|Hg(ω)|');title('(d)')
%============================================================================================
figure(2)   %图7.3.2
plot(wk/pi,Hgw);axis([0,1,-0.2,1.2])
hold on;
plot(wk1,Aw,':');plot(wk2,Ak,'o');
%=============================================================================================
%采样点数N=65
N=75;
Np=fix(wc/(2*pi/N));Ns=N-2*Np-1; %Np+1为通带[0,wc]上采样点数,Ns为阻带[wc,2*pi-wc]上采样点数
Ak=[ones(1,Np+1),zeros(1,Ns),ones(1,Np)];     %幅度采样向量A(k)
wk2=[0:N-1]*2/N;
plot(wk2,Ak,'.');
thetak=-pi*(N-1)*(0:N-1)/N;                 %相位采样向量θ(k) 
Hk=Ak.*exp(j*thetak);    %构造频域采样向量H(k)
hn=real(ifft(Hk));      %h(n)=IDFT[H(k)]
Hw=fft(hn,1024);
wk=2*pi*[0:1023]/1024;
Hgw=Hw.*exp(j*wk*(N-1)/2);
%=====================================================
%计算N=75时的通带最大衰减rp15和阻带最小衰减rs15
rp75=max(20*log10(abs(Hgw)))
hgmin=min(real(Hgw))
rs75=20*log10(abs(hgmin))
%========================================================
plot(wk/pi,Hgw);
xlabel('ω/π');ylabel('Hg(ω)')