经典谱估计

经典谱估计

周期图法

把随机信号$x(n)$的N点观察数据$x_N(n)$视为一个能量有限信号，直接取$x_N(n)$的傅里叶变换，得到$X_N(e^{jw})$，然后再取其幅值的平方，除以N，作为对$x(n)$真实的功率谱$P(e^{jw})$的估计。
$${\hat{P}}_{PER}(w)=\frac{1}{N}|X_N(w){|}^2$$

        [Pxx,f] = periodgram(xn,window,nfft,fs)
        window: rectwin(length(xn))     hamming(length(xn))  .....

自相关法

由$x_N(n)$估计出自相关函数$\hat{r}(m)$，然后对$\hat{r}(m)$求傅里叶变换得到$x_N(n)$的功率谱，记为${\hat{P}}_{BT}(w)$，并以此作为对$P(w)$的估计。
$${\hat{P}}_{BT}(w)=\sum _{m=-M}^M\hat{r}(m)e^{-jwm},|M|\le N-1$$

         cxn = xcorr(xn,'unbiased');
         P = fft(cxn);

Bertlett法

对L个具有相同均值$\mu$和方差${\sigma }^2$的独立随机变量$X_1,X_2,\cdots ,X_L$，新随机变量$X=(X_1+X_2+\cdots +X_L)/L$的均值也是$\mu$，但方差变为${\sigma }^2/L$，小了L倍。由此得到改善${\hat{P}}_{PER}(w)$方差特性的一个有效方法。将采样数据$x_N(n)$分为L段，每段的长度都是M，即N=LM，第i段数据加矩形窗后，变为$${\hat{P}}_{PER}^i(w)=\frac{1}{M}\sum _{i=1}^L{\left|\sum _{n=0}^{M-1}{x}_N^i(n)e^{-iwn}\right|}^2$$，把${\hat{P}}_{PER}^i(w)$对应相加，再取平均，得到平均周期图。
$${\bar{P}}_{PER}(w)=\frac{1}{L}\sum _{i=1}^L{\hat{P}}_{PER}^i(w)=\frac{1}{ML}\sum _{i=1}^L{\left|\sum _{n=0}^{M-1}{x}_N^i(n)e^{-iwn}\right|}^2$$

Welch法

Welch法是对Bartlett法的改进。在对$x_N(n)$分段时，可允许每一段数据有部分的交叠；每一段的数据窗口可以不是矩形窗口，例如使用汉宁窗或汉明窗，记为$d_2(n)$。这样可以改善由于矩形窗边瓣较大引起的频谱失真。
$${\bar{P}}_{PER}^i(w)=\frac{1}{MU}{\left|\sum _{n=0}^{M-1}{x}_N^i(n)d_2(n)e^{-jwn}\right|}^2$$
$$U=\frac{1}{M}\sum _{n=0}^{M-1}{d}_2^2(n)$$

          [Pxx,f] = pwelch(xn,window,noverlap,nfft,fs,range)

代码

%谱估计 2023.09.19
close all;clear;
N = 1024;
fs = 1000;
n = (0:N-1)/fs;
xn = sin(2* pi*100*n) + sin(2* pi*200*n) + 0.001*randn(size(n));

%周期图法
[Pxx,f] = periodogram(xn,hamming(N),1024,fs);
%相关法
cxn = xcorr(xn,'unbiased');
P = fftshift(fft(cxn));
%Welch
[Pxx2,f2] = pwelch(xn,400,200,N,fs,'onesided'); %每段400，有200重叠

figure;
subplot(1,3,1);plot(f,10*log10(abs(Pxx)));title('周期图')
subplot(1,3,2);plot((-N+1:N-1)/N*fs/2,10*log10(abs(P)));title('相关法')
subplot(1,3,3);plot(f2,10*log10(abs(Pxx2)));title('Welch')

时频分析

问题背景：在处理声致水面微动信号问题时，一般处理流程是提取水面振动信号$y(t)$以后，通过周期图等频谱估计方法估计声源频率。
当声源频率随时间变化时，我们不仅想要知道频率，更想知道频率随时间的变化，此时就需要用到时频分析。

spectrogram