[編者按]傳感器和信號(hào)處理僅一線之隔，信號(hào)的前后端合理搭配，是我們更準(zhǔn)確地感知這個(gè)世界的一種基本態(tài)度和方式。

FIR頻域重疊相加法

還記得我們（此處有重復(fù)之嫌）之前的發(fā)文《FIR連續(xù)采樣分段卷積時(shí)域重疊相加法》？不過(guò)那是在時(shí)域處理的模擬和仿真。這次我們的內(nèi)容是用C++在頻域?qū)崿F(xiàn)的濾波卷積法，仍然是重疊相加法，屆時(shí)大家可以比較一下兩種方式的差異。

基本是通過(guò)兩個(gè)處理過(guò)程。

（1）頻域的重疊相加法示意圖再拿來(lái)用一下。如下圖所示[1]。

（2）再借用一下的時(shí)域卷積經(jīng)傅里葉FFT變換后，在頻域成為對(duì)應(yīng)的相乘；然后再通過(guò)IFFT將中間結(jié)果轉(zhuǎn)換回時(shí)域時(shí)序結(jié)果。

讓我們直接跳進(jìn)話題，先看模擬測(cè)試結(jié)果，后看C++源碼。

模擬情節(jié)設(shè)定

50Hz選頻濾波，信號(hào)中混有110Hz和210H在的干擾信號(hào)和幅值為1的直流DC。

模擬信號(hào)及其頻譜的輸出請(qǐng)查看我們前面的文章。這里的代碼只提供將模擬信號(hào)進(jìn)行了頻域重疊相加處理，生成的濾波前后模擬信號(hào)和被濾波處理后的數(shù)據(jù)波形的比較（見(jiàn)下圖）。

還記得我們（此處重復(fù)）之前用C++來(lái)模擬時(shí)域處理的濾波模擬程序嗎？

你又猜對(duì)了，又是那個(gè)濾波器，又被用上了！但，是不同的實(shí)現(xiàn)處理方式。

濾波處理之前的波形和頻譜圖

濾波之后，直流和其他頻率的信號(hào)已經(jīng)不見(jiàn)，只留下50Hz的正弦波（見(jiàn)下圖）。

頻域重疊相加濾波前后的波形比較

圖由csv文件處理后生成。又見(jiàn)此圖，是不是有熟悉的感覺(jué)？

頻域連續(xù)濾波模擬和驗(yàn)證C++源碼

/*
Project Name: Demonstration & Validation for signal filtering via the 
"Overlap-Add" method implemented through FFT/IFFT based on Fast Fourier Transform (FFT)
based on the Cooley-Tukey algorithm.
Copyright(c)2024AmphenolSensors
Author: L.L.


This software is provided 'as-is', without any express or implied
warrantyandforsimulation/demostrationpurpose.Innoeventwilltheauthorsbeheldliableforanydamages
arising from the use of this software.


Permission is granted to anyone to use this software for any purpose,
including commercial applications, and to alter it and redistribute it
freely, subject to the following restrictions:


The origin of this software must not be misrepresented; you must not
claim that you wrote the original software. If you use this software
in a product, an acknowledgment in the product documentation would be
appreciated but is not required.
*/
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 


const double PI = 3.1415926535897932384;
const double SAMPLE_RATE = 1000.0;      // 1000Hz
const int N = 1024;                     // 假設(shè)采樣分段數(shù)為1024


using namespace std; // 聲明使用std命名空間


#define SEL_FFT     0
#define SEL_IFFT    1
#define SEL_FFTIFFT 2


#define SIM_CYCLE_CNT   3.0


// 模擬信號(hào)函數(shù)
vector> generateSignal(double sampleRate, double seconds)
{
    size_t signal_len = (size_t)(sampleRate * seconds);
    vector> signal(signal_len); // 定義模擬信號(hào)的數(shù)組長(zhǎng)度
    for (size_t i = 0; i < signal_len; ++i)
    {
        // 包含50Hz和110Hz,210Hz信號(hào),DC
        signal[i].real(1+ sin((2 * PI * (double)i * 50) / sampleRate) + sin((2 * PI * (double)i * 210) / sampleRate) + sin((2 * PI * (double)i * 110) / sampleRate));
        signal[i].imag(0);
    }
    return signal;
}


// 基于Cooley-Tukey算法的FFT
void fft2(vector> &x)
{
    size_t n = x.size();
    if (n <= 1)
        return;


    // 把序列分為奇偶兩部分
    vector> even(n / 2), odd(n / 2);
    for (size_t i = 0; 2 * i < n; i++)
    {
        even[i] = x[i * 2];
        odd[i] = x[i * 2 + 1];
    }


    // 遞歸解決每一個(gè)序列
    fft2(even);
    fft2(odd);


    // 結(jié)合步驟
    for (size_t i = 0; 2 * i < n; i++)
    {
        complex t = polar(1.0, -2 * PI * (double)i / (double)n) * odd[i];
        x[i] = even[i] + t;
        x[i + n / 2] = even[i] - t;
    }
}
// 基于Cooley-Tukey算法的IFFT
void ifft2(vector> &x)
{
    size_t n = x.size();
    if (n <= 1)
        return;


    vector> even(n / 2), odd(n / 2);
    for (size_t i = 0; 2 * i < n; i++)
    {
        even[i] = x[i * 2];
        odd[i] = x[i * 2 + 1];
    }


    ifft2(even);
    ifft2(odd);


    for (size_t i = 0; 2 * i < n; i++)
    {
        complex t = polar(1.0, 2 * PI * (double)i /(double) n) * odd[i];
        x[i] = even[i] + t;
        x[i + n / 2] = even[i] - t;
    }
    // 最后除以n
    if (n > 1)
    {
        for (size_t i = 0; i < n; i++)
        {
            x[i] /= 2;
        }
    }
}


// 將數(shù)據(jù)寫入文件
void writeToFile(const vector> &Original_signal, const vector> &Proceeded_Signal, const string &filename, int Type_sel)
{
    ofstream file(filename);
    if (Type_sel==SEL_FFTIFFT)
    {
        file << "time(s)" << ", " << "Original Signal"<< ", " << "Filtered Signal" << "
";
        for (size_t i = 0; i < Original_signal.size(); i++)
        {
            file << (double)i / SAMPLE_RATE << ", " < b{0.0010175493084400998, 0.0010954624020866333, 0.001080635650435545, 0.0009293052645812359,
                   0.0005868808563577278, -8.138309855847798e-19, -0.0008644147524968251, -0.0019966389877814107,
                   -0.003323586744207458, -0.004696461345361978, -0.005892320462621699, -0.006633249964255378,
                   -0.006623614506478284, -0.005601944833604465, -0.0034001773970723163, -7.334366341273803e-18,
                   0.004425290874832446, 0.00949426225087417, 0.014634010415364655, 0.019132982942933127,
                   0.022226796444847933, 0.023207550009729024, 0.021541722692400025, 0.01697833945185371,
                   0.009628503914736117, -6.755395515820625e-18, -0.01102370844120733, -0.02226281209657117,
                   -0.032372473621654914, -0.04001099412924139, -0.04402269970024527, -0.043609484958132556,
                   -0.03846490807520255, -0.028848803480728435, -0.015588116829396594, -9.10410551538968e-18,
                   0.016255406162706088, 0.031374390998733945, 0.04363491329762711, 0.051616779739690075,
                   0.05438594145724075, 0.051616779739690075, 0.04363491329762711, 0.031374390998733945,
                   0.016255406162706088, -9.10410551538968e-18, -0.015588116829396594, -0.028848803480728435,
                   -0.03846490807520255, -0.043609484958132556, -0.04402269970024527, -0.0400109941292414,
                   -0.032372473621654914, -0.022262812096571168, -0.01102370844120733, -6.755395515820627e-18,
                   0.009628503914736117, 0.016978339451853702, 0.021541722692400025, 0.023207550009729034,
                   0.022226796444847933, 0.01913298294293312, 0.014634010415364655, 0.009494262250874175,
                   0.004425290874832446, -7.3343663412738e-18, -0.0034001773970723163, -0.005601944833604469,
                   -0.006623614506478284, -0.006633249964255374, -0.005892320462621699, -0.00469646134536198,
                   -0.003323586744207458, -0.001996638987781409, -0.0008644147524968251, -8.138309855847805e-19,
                   0.0005868808563577278, 0.0009293052645812359, 0.001080635650435545, 0.0010954624020866333,
                   0.0010175493084400998};


    // (1) Resize filter coefficients
    vector> H(b.size());


    for(size_t i=0; i< b.size(); i++)
    {
        H[i] = complex(b[i],0);
    }
    H.resize(N, complex(0.0, 0.0));


    // (2)Generate simulation data sequences
    size_t DataSeg_len_L = N - b.size() + 1;                // Data segmeng Length = L


    double daq_duration = (double)DataSeg_len_L * SIM_CYCLE_CNT / SAMPLE_RATE;
    vector>Original_signal = generateSignal(SAMPLE_RATE, daq_duration);     // Generate data sequence, L * 3


    // (3-1) Define a 2-D vector，3 columns, to simulate DAQ and filtering process for 3 rounds
    vector>> seg_Dates(3); 
    // (3-2) Initialize data segment
    for (size_t i = 0; i < seg_Dates.size(); i++)
    {
        seg_Dates[i].resize(DataSeg_len_L,complex(0.0, 0.0));
        // devide Original_signal into 3 parts to simulate consequent DAQ for 3 cycles
        for(size_t j=0; j(0.0, 0.0));                     // resize each data segment to N = L + M - 1
    }


    // (4) Start to FFT/IFFT and generate involution result
    vector> Filtered_signal(DataSeg_len_L * (size_t)(SIM_CYCLE_CNT) +  b.size() -1 );     //L * 3 + (M - 1)
    fft2(H);  


//(4-1)Simulate3cyclesofinvolution:L*3
    for(size_t i=0; i< seg_Dates.size(); i++)
    {
        fft2(seg_Dates[i]);
        for(size_t j = 0; j< seg_Dates[i].size(); j++)
        {
            seg_Dates[i][j] = seg_Dates[i][j] * H[j];
        }


        ifft2(seg_Dates[i]);
        for(size_t k = 0; k< seg_Dates[i].size(); k++)
        {
            Filtered_signal[i*DataSeg_len_L + k] += seg_Dates[i][k];
        }
    }


????//?(5)?Save?Origianl?signal & result?(data?after?filtering)?into?csv?file
    writeToFile(Original_signal, Filtered_signal,"output_Filtered2.csv", SEL_FFTIFFT);


    return 0;
}

時(shí)間倉(cāng)促，有些功能和細(xì)節(jié)并沒(méi)有考慮太多，這里功能驗(yàn)證是第一。

FFT/IFFT是基于庫(kù)里-圖基的算法，各位可以選用其他的來(lái)優(yōu)化替代；

有些參數(shù)可以單獨(dú)變，有些參數(shù)卻是關(guān)聯(lián)的。