hc741b.m

书籍中matlab第二部分代码

%《数字信号处理教程——MATLAB释义与实现》第七章例7.4.1b程序hc741b
% 取不同样本数时的幅频特性
% 电子工业出版社出版  陈怀琛编著 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] = myfreqz(h,[1]);    % 检验设计出的滤波器的频率响应
[Ar,ww,type,L0]=amplres(h);             % 检验设计出的滤波器的符幅特性
dw = 2*pi/1000;
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]; % 绘制理想幅特性的频率和幅度数据
plot(ww,mag,wl(1:N),A(1:N),'.');            % 画图
set(gcf,'color','w'),hold on
axis([0,4,0,1.2]); title('幅频特性')
xlabel('频率（单位：弧度/秒）'); ylabel('|H|')
set(gcf,'color','w');                % 置图形背景色为白