|
下面的matlab程序分别使用周期图法、相关函数法以及AR谱方法计算信号的功率谱 - % power spectrum estimated
- clear all;
- clc;
- close all;
- Fs=1000; % 采样频率
- nfft = 1024; % fft计算点数
- %产生含有噪声的序列
- n=0:1/Fs:1;
- xn=cos(2*pi*100*n)+3*cos(2*pi*200*n)+(2*randn(size(n)));
- subplot(2,2,1);plot(xn);title('加噪信号');grid on
- % 周期图法
- window=boxcar(length(xn)); %矩形窗
- [psd1,f]=periodogram(xn,window,nfft,Fs); %直接法
- psd1 = psd1 / max(psd1);
- subplot(2,2,2);plot(f,10*log10(psd1+0.000001));title('周期图法');grid on
- % 自相关法
- cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
- CXk=fft(cxn,nfft);
- psd2=abs(CXk);
- index=0:round(nfft/2-1);
- k=index*Fs/nfft;
- psd2 = psd2 / max(psd2);
- psd2=10*log10(psd2(index+1)+0.000001);
- subplot(2,2,3);plot(k,psd2);title('自相关法');grid on
- % AR谱
- psd3 = pyulear(xn, Fs, nfft);
- psd3=psd3/max(psd3);
- psd3 = psd3 / max(psd3);
- index=0:round(nfft/2-1);
- psd3=10*log10(psd3(index+1)+0.000001);
- subplot(2,2,4);plot(k, psd3);title('AR谱估计');grid on;
| 
楼主
