平稳序列自相关和功率谱密度的估计.docVIP

实平稳序列自相关和功率谱密度的估计 1.直接法：带偏的一致估计 ① 其中N是序列的长度。 2.间接法：利用DTFT的卷积性质 序列{x(n)}{y(n)}长度分别为L1,L2，则对上式求IDTFT便得到x(n),y(n)线性卷积的结果序列，长度为M =L1+ L2 -1。 由于DTFT是周期连续的，可以由DFT近似表示其在一个周期内的离散频谱。即 因为在DFT中，频域上M点乘积对应时间域上两个长为M的序列的循环（周期）卷积。 ＝ （表示M点循环卷积） 当M =L1+ L2 -1时，循环卷积在一个周期内的主值就是长分别为L1和L2的序列（分别在尾部加L2-1和L1-1个零）的线性卷积。 所以，在用间接法求自相关时，M=2N-1。利用DFT的性质x(-n)→X(-k)，以及x(n)是实序列这一事实，利用频谱的Hermitian性质X((-k)=X*(k)，得 , 时间域上，在序列x(n)的尾部加N-1个0，新序列记作x1(n)。对x1(n)求DFT（使用FFT算法），得到2N-1点的X(k)。将其取模平方后求IDFT（使用FFT算法求ifft），便得到了循环卷积x1(-n) x1(n)在一个周期（M个点）内的值。它等于x(-n)和x(n)的线性卷积。在此基础上除以序列长度N（注意：不是除以频域点数，也不是除以M），便得到自相关函数R(n)。 Matlab代码(间接法求自相关)： xn=randn(1,1000); N=length(xn); Xk=fft(xn, 2*N-1); %M=2*N-1,其实M也可以取更大的值，以方便快速运算DFT Rn=ifft(abs(Xk).^2)/N; Rn=[Rn(end-N+2:end) Rn(1:N)]; %这种做法适用于fft点数多于所需的2*N-1个的时候 n=(-N+1:N-1); figure;plot(n, Rn); 图1 间接法求0均值单位方差高斯噪声过程的自相关函数 经典谱估计也分为直接法和间接法。直接法利用序列的DFT模平方除以频域序列点数M（这个点数可以和时域信号点数N不同）；间接法利用自相关函数进行DFT。DFT都采用FFT算法。用matlab中的pwelch要乘以fs(估计双边谱)，或fs/2(单边谱)。 经典功率谱估计的Matlab例程： %class_psd.m %利用经典方法估计PSD， %直接法直接从序列的DFT取模平方之后除以DFT的点数来估计PSD。 %间接法先计算序列的自相关函数，然后作DFT，作为PSD的估计。 %welch方法是修正的周期图法（直接法），要加窗，重叠，叠加 N = 1024; t=linspace(0,1,N); signal=4*sin(2*pi*50*t)+5*sin(2*pi*200*t); %signal = 0; Ts = 1e-3; fs = 1/Ts; noise=randn(size(t)); symbol=signal+noise; Y=fft(symbol,N*2); f=(0:N)*fs/(2*N); P1=Y.*conj(Y)/(2*N); %直接法，要除的是FFT的点数 subplot(311); plot(f,10*log10(P1(1:N+1))); xlabel(frequency) ylabel(power(dB)) title(直接法,序列DFT模平方再除以DFT点数) %间接法 cxn=xcorr(symbol,unbiased); %计算序列的自相关函数，2N-1点，中心点对应R(0)。 P2=abs(fft(cxn, 2*N)); %采用2N个点的FFT subplot(312); f=(0:N-1) * fs/(2*N); %fs/(2N)是频率分辨率 plot(f,10*log10(P2(1:N))); xlabel(frequency(Hz)) ylabel(power(dB)) title(间接法，利用自相关的DFT); %加窗法 window2=blackman(100); %blackman窗 noverlap=32; %数据重叠 range=onesided; %频率间隔为[0 fs/2]，只计算一半的频率 [P3(1:N/2+1),f]=pwelch(symbol,window2,noverlap,N,fs,range); P3=P3.*fs/2; %单边，要乘以fs/2 plot_P3=10*log10(P3(1:N/2+1)); subplo