從多個方面詳解pwelch函數

一、pwelch函數的用法

pwelch函數是一種信號處理函數，可以用於估計與輸入信號相關的功率譜密度。該函數主要用於雜訊分析、信噪比(SNR)分析以及頻域特徵提取等領域。

使用pwelch函數時，需要輸入一組時域信號，輸出的是功率譜密度估計。

Matlab函數原型：
[Pxx, F] = pwelch(x, window, noverlap, NFFT, Fs);
C++實現：
void pwelch(const std::vector& x, std::vector& pxx, std::vector& f,
            int window, int noverlap, int nfft, double fs);

其中，x是輸入的時域信號；window是指定的窗函數類型及其長度；noverlap是指定重疊樣本數；NFFT是離散傅里葉變換(DFT)點數；Fs是採樣頻率。在C++實現中，輸出參數pxx和f分別是功率譜密度估計和對應的頻率坐標軸。

二、pwelch函數的讀音

pwelch函數的讀音為「皮威爾奇」。

三、pwelch函數的實現原理

pwelch函數的實現基於傅里葉變換和窗函數。窗函數有多種類型可選，比如矩形窗、漢寧窗、漢明窗等，它們的作用是在時域上減小信號的尾部效應，使頻域上的信號更加平滑。一般來說，窗口的長度越大，解析度越高，但能量泄漏的越大，頻譜圖會更模糊。而重疊樣本數決定著每一個時間窗之間的數據重合程度，增加重疊可以從一定程度上提高頻譜估計的準確性。

具體的，pwelch函數將輸入信號分為多個段，在每個段內，用窗函數挑選出所需要的長度並加以平方，計算每個段的傅里葉變換，並對所有段的頻譜進行平均處理，得到最終的功率譜密度估計。

四、pwelch函數的C++實現

以下是採用C++實現的pwelch函數代碼示例：

void pwelch(const std::vector& x, std::vector& pxx, std::vector& f,
            int window, int noverlap, int nfft, double fs) {
    // 計算窗口長度
    int nw = window;
    if (nw == 0) { // 如果未輸入窗口長度，則自動選擇
        double duration = x.size() / fs;
        nw = std::round(std::max(256.0, 2 * fs * duration - noverlap) / (1 + noverlap));
        nw = (nw % 2 == 0 ? nw : nw + 1); // 窗口長度需要為偶數
    }
    int nfftActual = nfft != 0 ? nfft : nw; // 計算實際的FFT點數

    // 初始化
    pxx.clear();
    f.clear();
    pxx.resize(nfftActual / 2 + 1);
    f.resize(nfftActual / 2 + 1);
    std::fill(pxx.begin(), pxx.end(), 0);
    std::fill(f.begin(), f.end(), 0);

    // 計算窗口函數
    std::vector windowFunc(nw);
    calculateWindowFunction(nw, windowFunc);

    // 計算分段數量
    int nseg = 1 + std::floor((x.size() - nw) / (nw - noverlap));

    // 計算每一段的功率譜密度估計
    std::vector segPower(nfftActual);
    for (int i = 0; i < nseg; i++) {
        // 計算起點和終點
        int start = i * (nw - noverlap);
        int end = start + nw;

        // 取出該段數據，並作窗函數處理
        std::vector segData(nw);
        for (int j = start, k = 0; j < end; j++, k++) {
            segData[k] = x[j] * windowFunc[k];
        }

        // 計算FFT，並提取出幅度譜的一半
        std::vector fftResult(2 * nfftActual);
        FFT::computeRealFFT(segData, fftResult);
        for (int j = 0; j <= nfftActual / 2; j++) {
            double re = fftResult[2*j];
            double im = fftResult[2*j+1];
            segPower[j] = 20.0 * log10(std::sqrt(re*re + im*im) / nw);
        }

        // 疊加各段的功率譜密度估計
        for (int j = 0; j <= nfftActual / 2; j++) {
            pxx[j] += segPower[j];
        }
    }

    // 對每一段的平均功率譜密度估計除以段數，然後取出對應的頻率坐標軸
    for (int i = 0; i <= nfftActual / 2; i++) {
        pxx[i] /= nseg;
        f[i] = i * fs / nfftActual;
    }
}

五、pwelch函數在Matlab中的應用

在Matlab中，pwelch函數的用法是相似的，比如下面的代碼：

% 產生隨機雜訊
Fs = 1000;
t = 0:(1/Fs):1;
x = rand(size(t));

% 處理數據並繪製功率譜密度圖
[Pxx, F] = pwelch(x, 1000, 500, 10000, Fs);
plot(F, Pxx);
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dB/Hz)');

六、pwelch函數與FFT的區別

FFT是一種快速傅里葉變換，常用於計算大量數據的DFT。FFT只能計算離散時間信號的傅里葉變換，輸出的是信號在頻域上的各頻率分量的幅度和相位信息。

pwelch函數是在FFT的基礎上，利用了窗函數等技術對輸入信號進行分段處理，並對多個段的頻譜估計值取均值得到更加平滑的功率譜密度估計。所以，pwelch函數在雜訊分析、SNR分析以及頻域特徵提取等領域具有很高的實用價值。