hc741.m

这是数字信号处理——matlab释义与实现一书的源代码

%《数字信号处理教程——MATLAB释义与实现》第七章例7.4.1程序hc741
% 用频率样本法设计FIR滤波器
% 电子工业出版社出版  陈怀琛编著 2004年9月
%
N=input('N= (N必须为奇数）');wc=input('wc=');     % 给出原始数据
N=N+mod(N+1,2);
N1=fix(wc/(2*pi/N));N2=N-2*N1-1;
A=[ones(1,N1+1),zeros(1,N2),ones(1,N1)];	% 符幅特性样本序列
theta = - pi*[0:N-1]*(N-1)/N;			% 相位特性样本序列
H=A.*exp(j*theta);				% 频率特性样本序列
h=real(ifft(H));				% 反变换求出脉冲序列，去掉运算误差造成的虚部
wp1=2*pi/N*fix(wc/(2*pi/N));ws1=wp1+2*pi/N;
[db,mag,pha,grd,w] = freqz_m(h,[1]);    % 检验设计出的滤波器的频率响应
[Ar,ww,type,L0]=amplres(h);            % 检验设计出的滤波器的符幅特性
dw = 2*pi/1000;                    % myfreqz函数中的频率分辨率
Rp = -min(db(1:fix(wp1/dw)+1))      % 实际的通带波动
As = -round(max(db(fix(ws1/dw)+1:501)))  % 最小阻带衰减
l = 0:N-1; wl = (2*pi/N)*l;        % 由频率样本下标换成频率样本值 
wdl=[0,wc,wc,2*pi-wc,2*pi-wc,2*pi]/pi;Adl=[1,1,0,0,1,1];    % 绘制理想幅特性的频率和幅度数据
subplot(2,2,1);plot(wl(1:N)/pi,A(1:N),'.',wdl,Adl);         % 绘图
axis([0,1,-0.1,1.1]); title('频率样本')
xlabel(''); ylabel('A(k)')
set(gca,'XTickMode','manual','XTick',chop([0,2/N,4/N,1],2))
set(gca,'YTickMode','manual','YTick',[0,1]); grid
subplot(2,2,2); stem(l,h,'.'); %axis([-1,N,-0.1,0.5])
title('脉冲响应');ylabel('h(n)');text(N+1,-0.1,'n')
subplot(2,2,3); plot(ww/pi,Ar,wl(1:N)/pi,A(1:N),'.');
axis([0,1,-0.2,1.2]); title('符幅响应')
xlabel('频率（单位：pi）'); ylabel('Ar(w)')
set(gca,'XTickMode','manual','XTick',chop([0,2/N,4/N,1],2))
set(gca,'YTickMode','manual','YTick',[0,1]); grid
subplot(2,2,4);plot(w/pi,db); axis([0,1,-50,10]); grid
title('幅度响应'); xlabel('频率（单位：pi）');
ylabel('分贝数');
set(gca,'XTickMode','manual','XTick',chop([0,2/N,4/N,1],2))
set(gca,'YTickMode','Manual','YTick',[-As;0]);
set(gca,'YTickLabelMode','manual','YTickLabels',['As';' 0'])
set(gcf,'color','w');                % 置图形背景色为白