matlab 数字信号处理 周期图 程序.docVIP

1、假设一平稳随机信号为，其中是均值为0，方差为1的白噪声，数据长度为1024。 （1）、产生符合要求的和； （2）、给出信号x(n)的理想功率谱； （3）、编写周期图谱估计函数，估计数据长度N=1024及256时信号功率谱，分析估计效果。 （4）、编写Bartlett平均周期图函数，估计当数据长度N=1024及256时，分段数L分别为2和8时信号的功率谱，分析估计效果。 解:思路 在matlab中提供的有 randn（m.n）函数，其为均值为零，方差为1的函数，所以w(n)可以通过随机序列randn(1,N)来产生，x(n)可以通过对 w(n)滤波产生，也可以直接由递推式迭代产生。 由于线性系统的输出功率谱等于输入功率谱乘以传递函数模的平方，X(n)可以看做w(n) 通过一线性系统的输出， H(z)=1/(1-0.8z) 所以x(n)的理想功率谱 直接产生 Matlab 程序： clear;close all; F=500; %采样率 N=1024; %观测数据 subplot(2,1,1); w=sqrt(1)+randn(1,N); plot(w); xlabel（观察次数）；ylabel(功率db)；title(白噪声w的分布情况); subplot(2,1,2); x=[w(1) zeros(1,N-1)]; %初始化x(n),长度1024,x(1)=w(1) for i=2:N x(i)=0.8*x(i-1)+w(i); %迭代产生观测数据x(n) end plot(x); xlabel(观测次数);ylabel(x 的分布情况);title(x的分布情况); %% 理想功率谱 cn=xcorr(x,x);%自相关函数 Cn=fft(cn,N);%快速傅里叶变化 Pxx=abs(Cn);%取绝对值 index=0:round(N/2-1); k=index*Fs/N; pxx=10*log10(Pxx(index+1));%转换成db figure; plot(k,pxx); xlabel(频率);ylabel(功率db); title(理想功率谱)； %%周期图谱 %1024 Pxx1=abs(fft(x)).^2/N; pxx1=10*log10(Pxx1(index+1)); figure; plot(k,pxx1); title(周期图1024个点); xlabel(频率);ylabel(功率db)； % 周期图256个观测点 x1=x(1:4:N); Pxx2=abs(fft(x1,1024)).^2/N; pxx2=10*log10(Pxx2(index+1)); figure; plot(k,pxx2); title(周期图256)；xlabel(频率); ylabel(功率db); L=2; x1=x(1:L:N); x2=x(2:L:N); pxx2_1=abs(fft(x1,1024)).^2/length(x1); pxx2_2=abs(fft(x2,1024)).^2/length(x2); pxx_2=(pxx2_1+pxx2_2)/L; figure; subplot(2,1,1); plot(k,10*log10(pxx_2(index+1))); title(N=1024,L=2时的周期图); xlabel(频率); ylabel(功率db); L=2; x1=x(1:L:N); x2=x(2:L:N); pxx2_1=abs(fft(x1,1024)).^2/length(x1); pxx2_2=abs(fft(x2,1024)).^2/length(x2); pxx_2=(pxx2_1+pxx2_2)/L; figure; subplot(2,1,1); plot(k,10*log10(pxx_2(index+1))); title(N=256,L=2时的周期图); xlabel(频率); ylabel(功率db); 从上图可以看出， 周期图法得到的功率谱估计， 谱线的起伏较大， 即估计所得的均方误差较大。当 N加时，摆动的频率加快，而摆动的幅度变化不大。且 N=1024时，谱的分辨率较N=256时大。 采用平均处理后， 谱线上下摆动的幅度减小（即均方误差有所降低），曲线的平滑性也较周期图法好。相同时，分段越多，方差越小，曲线越平滑。这是因为， N一定时L加大，每一段的数据量就会相应减少，因此估计方差减小，偏移加大，从而分辨率降低。 N一定时，L 与每段数据量相互矛盾， 需择中选取。 综上， Bartlet