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

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);%转换成dbfigure;plot(k,pxx); xlabel(频率);ylabel(功率db); title(理想功率谱)；%周期图谱%1024Pxx1=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 N=1024 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); L1=8;%当数据长度为8时x3=zeros(L1,N/L1);for i=1:L1x3(i,:)=x(i:L1:N);end pxx3=zeros(L1,N);for i=1:L1pxx3(i,:)=abs(fft(x3(i,:),1024).2/length(x3(i,:);end for i=1:1024Pxx3(i)=sum(pxx3(:,i)/L1;end subplot(2,1,2);plot(k,10*log10(Pxx3(index+1);title(N=1024,L=8时的周期图); xlabel(频率);ylabel(功率db);当长度为256时，指针函数发生变化clear;close all;Fs=500; %采样率N=256; %观测数据w=sqrt(1)+randn(1,N); %x=w(1) zeros(1,N-1);for i=2:Nx(i)=0.8*x(i-1)+w(i);endindex=0:round(N/2-1);k=index*Fs/N; 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); L1=8;%当数据长度为8时x3=zeros(L1,N/L1);for i=1:L1x3(i,:)=x(i:L1:N);end pxx3=zeros(L1,N);for i=1:L1pxx3(i,:)=abs(fft(x3(i,:),1024).2/length(x3(i,:);end for i=1:1024Pxx3(i)=sum(pxx3(:,i)/L1;end subplot(2,1,2);plot(k,10*log10(Pxx3(index+1);title(N=256,L=8时的周期图); xlabel(频率);ylabel(功率db);结果:理想功率谱：长度为256跟长度为1024时候的周期图谱：N=1024时，当L=2与L=8时 的周期图当N的长度为256时，L=2与L=8时的周期图谱 从上面的图像可以看出， 周期图法得到的功率谱估计， 谱线的起伏较大， 即估计所得的均方误差较大。当 N加时，摆动的频率加快，而摆动的幅度变化不大。且 N=1024时，谱的分辨率较N=256时大。 采用平均处理后， 谱线上下摆动的幅度减小（即均方误差有所降低），曲线的平滑性也较周期图法好。N相同时，分段越多，方差越小，曲线越平滑。这是因为， N一定时L加大，每一段的数据量就会相应减少，因此估计方差减小，偏移加大，从而分辨率降低。 N一定时，L 与每段数据量相互矛盾， 需择中选取。 综上， Bartlett 法相对于周期图法来说，较好地减小了估计误差。 2、假设均值为0，方差为1的白噪声中混有两个正弦信号，该正弦信号的频率分别为100Hz和110Hz，信噪比分别为10dB和30dB，初始相位都为0，采样频率为1000Hz。（1）、采用自相关法、Burg法、协方差法、修正协方差法估计功率谱，分析数据长度和模型阶次对估计结果的影响（可采用MATLAB自带的功率谱分析函数）。（2）、调整正弦信号信噪比，分析信噪比的降低对估计效果的影响。Matlab源程序clc;clear;close all;fs=1000;%采样率1000N=1;%改变数据长度p=50;%AR模型阶数nfft=512;%fft长度t=0:1/fs:N;wn=sqrt(1)+randn(1,N*fs+1); %白噪声,均值0,方差1s1=sqrt(20)*sin(2*pi*100*t); %正弦信号1,信噪比10dbs2=sqrt(2000)*sin(2*pi*110*t); %正弦信号2,信噪比30dbx=s1+s2+wn; %观测数据%figure,plot(t,x);x1=xcorr(x,biased);Pxx,f=pyulear(x1,p,nfft,fs); %Yule-Walker方程figure,plot(f,10*log10(Pxx);grid on;title(自相关法);Pxx1,f1=pcov(x,p,nfft,fs);figure,plot(f1,10*log10(Pxx1);grid on;title(协方差法);Pxx2,f2=pmcov(x,p,nfft,fs);figure,plot(f2,10*log10(Pxx2);grid on;title(修正协方差法);Pxx3,f3=pburg(x,p,nfft,fs); figure,plot(f3,10*log10(Pxx3);grid on;title(Burg法);结果：长度为N=1，p=50时的图像，长度N=5，p=50是的图片N=1， p=100时的图谱N=1，p=80 的图谱分析数据长度和模型阶次对估计结果的影响：从上面的图像可以看出，当数据的长度