hc434a.m

某本matlab书籍中的代码

%《数字信号处理教程——MATLAB释义与实现》第四章例4.3.4a程序hc434a
% 计算无限长序列频谱时截取长度的自动选择程序
% 电子工业出版社出版  陈怀琛编著 2004年9月
%
T=0.5; a=1; b=100; beta=1;      % 给定初始数据
while b>beta                    % 判断是否应结束循环运算
    N1=2^a;n1=0:N1-1;           % 确定数据长度N1
    x1=0.7.^n1; X1=fft(x1);     % 求长度N1的序列x及其FFT X1   
    N2=2*N1; n2=0:N2-1;         % 数据长度加倍为N2=2*N1
    x2=0.7.^n2;X2=fft(x2);      % 求长度N2的序列x2及其FFT X2
    k1p=0:N1/2;k2p=2*k1p;       % 确定两序列对应点的下标
    d=max(abs(X1(k1p+1)-X2(k2p+1))); % 求对应点上FFT的误差
    Xm=max(abs(X1(k1p+1)));     % 求X1幅特性的最大值
    b=d/Xm*100;                 % 求相对误差的百分数
    a=a+1;                      % 序列加长一倍
N2,b                          % 结束循环后显示达到要求的长度N2和相对误差b
k=floor(-N2/2+0.5:N2/2-0.5);dw=2*pi/(N2*T);     % 绘图用的乃奎斯特频率范围
subplot(1,2,1) 
plot(k*dw,abs(fftshift(X2))),title('(a)振幅频谱'),hold on % 画幅特性
xlabel('频率'),ylabel('幅度')
set(gcf,'color','w')              % 置图形背景色为白
subplot(1,2,2) 
plot(k*dw,angle(fftshift(X2))),title('(b)相位频谱'),hold on % 画相特性
xlabel('频率'),ylabel('相位')
line([-10,10],[0,0])
end