ep731.m

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

%《数字信号处理--原理、实现及应用》第七章 例7.3.2 程序：ep732.m
% 电子工业出版社出版  高西全  丁玉美  阔永红合著 2006年8月
% 用频率采样法设计FIR低通滤波器
T=input('T= ')  %键入过渡采样值T
datB=pi/16;wc=pi/3;   %过渡带宽度pi/16,通带截止频率为pi/3；
m=1;
N=(m+1)*2*pi/datB+1;     %按（7.3.15）式估算采样点数N
N=N+mod(N+1,2);          %确保h(n)长度N为奇数
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)];     %N为奇数，幅度采样向量A(k)=A(N-k)
Ak(Np+2)=T;Ak(N-Np)=T;      %加一个过渡采样
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);    %计算频率响应函数:DFT[h(n)]
wk=2*pi*[0:1023]/1024;
Hgw=Hw.*exp(j*wk*(N-1)/2);  %计算幅度响应函数Hg(ω)
%=====================================================
%计算通带最大衰减Rp和阻带最小衰减Rs
Rp=max(20*log10(abs(Hgw)))
hgmin=min(real(Hgw));Rs=20*log10(abs(hgmin))
%==============================================================================================
%以下绘图7.3.4
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('Hg(k)');
hold on
wk2=[0:N-1]*2/N;plot(wk2,Ak,'.');
n=0:N-1;
subplot(3,2,3);stem(n,hn,'.');axis([0,65,-0.2,0.5]);xlabel('n');ylabel('h(n)');
subplot(3,2,5);plot(wk/pi,20*log10(abs(Hgw)));axis([0,1,-60,3]);grid on
xlabel('ω/π');ylabel('20lg|Hg(ω)|');