资源描述:
《matlab 数字信号处理 周期图 程序》由会员上传分享,免费在线阅读,更多相关内容在行业资料-天天文库。
1、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)滤波产生,也可以直接由
2、递推式迭代产生。由于线性系统的输出功率谱等于输入功率谱乘以传递函数模的平方,X(n)可以看做w(n)通过一线性系统的输出,H(z)=1/(1-0.8z)所以x(n)的理想功率谱直接产生Matlab程序:>>clear;closeall;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)];
3、%初始化x(n),长度1024,x(1)=w(1)fori=2:Nx(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
4、);>>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));>>fig
5、ure;>>plot(k,pxx2);>>title('周期图256');xlabel('频率');ylabel('功率db');%%L=2N=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)));
6、>>title('N=1024,L=2时的周期图');>>xlabel('频率');>>ylabel('功率db');>>L1=8;%当数据长度为8时x3=zeros(L1,N/L1);fori=1:L1x3(i,:)=x(i:L1:N);end>>pxx3=zeros(L1,N);fori=1:L1pxx3(i,:)=abs(fft(x3(i,:),1024)).^2/length(x3(i,:));end>>fori=1:1024Pxx3(i)=sum(pxx3(:,i))/L1;end>>subplot(2,1,2);plot(k,10*log1
7、0(Pxx3(index+1)));title('N=1024,L=8时的周期图');xlabel('频率');ylabel('功率db');当长度为256时,指针函数发生变化clear;closeall;Fs=500;%采样率N=256;%观测数据w=sqrt(1)+randn(1,N);%x=[w(1)zeros(1,N-1)];fori=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=ab
8、s(fft(x1,1024)).^2/length(x1);>>pxx2_2=abs(fft(