数字滤波器(6)—FIR频域连续滤波“重叠相加法”C++源码

[编者按]传感器和信号处理仅一线之隔，信号的前后端合理搭配，是我们更准确地感知这个世界的一种基本态度和方式。

FIR频域重叠相加法

还记得我们（此处有重复之嫌）之前的发文《FIR连续采样分段卷积时域重叠相加法》？不过那是在时域处理的模拟和仿真。这次我们的内容是用C++在频域实现的滤波卷积法，仍然是重叠相加法，届时大家可以比较一下两种方式的差异。

基本是通过两个处理过程。

（1）频域的重叠相加法示意图再拿来用一下。如下图所示[1]。

（2）再借用一下的时域卷积经傅里叶FFT变换后，在频域成为对应的相乘；然后再通过IFFT将中间结果转换回时域时序结果。

让我们直接跳进话题，先看模拟测试结果，后看C++源码。

模拟情节设定

50Hz选频滤波，信号中混有110Hz和210H在的干扰信号和幅值为1的直流DC。

模拟信号及其频谱的输出请查看我们前面的文章。这里的代码只提供将模拟信号进行了频域重叠相加处理，生成的滤波前后模拟信号和被滤波处理后的数据波形的比较（见下图）。

还记得我们（此处重复）之前用C++来模拟时域处理的滤波模拟程序吗？

你又猜对了，又是那个滤波器，又被用上了！但，是不同的实现处理方式。

滤波处理之前的波形和频谱图

滤波之后，直流和其他频率的信号已经不见，只留下50Hz的正弦波（见下图）。

频域重叠相加滤波前后的波形比较

图由csv文件处理后生成。又见此图，是不是有熟悉的感觉？

频域连续滤波模拟和验证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) 2024 Amphenol Sensors
Author: L.L.

This software is provided 'as-is', without any express or implied
warranty and for simulation/demostration purpose. In no event will the authors be held liable for any damages
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 <complex>
#include <cmath>
#include <stdio.h>
#include <iostream>
#include <vector>
#include <fstream>
#include <math.h>
#include <algorithm>

const double PI = 3.1415926535897932384;
const double SAMPLE_RATE = 1000.0;      // 1000Hz
const int N = 1024;                     // 假设采样分段数为1024

using namespace std; // 声明使用std命名空间

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

#define SIM_CYCLE_CNT   3.0

// 模拟信号函数
vector<complex<double>> generateSignal(double sampleRate, double seconds)
{
    size_t signal_len = (size_t)(sampleRate * seconds);
    vector<complex<double>> signal(signal_len); // 定义模拟信号的数组长度
    for (size_t i = 0; i < signal_len; ++i)
    {
        // 包含50Hz和110Hz,210Hz信号,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<complex<double>> &x)
{
    size_t n = x.size();
    if (n <= 1)
        return;

    // 把序列分为奇偶两部分
    vector<complex<double>> 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];
    }

    // 递归解决每一个序列
    fft2(even);
    fft2(odd);

    // 结合步骤
    for (size_t i = 0; 2 * i < n; i++)
    {
        complex<double> 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<complex<double>> &x)
{
    size_t n = x.size();
    if (n <= 1)
        return;

    vector<complex<double>> 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<double> 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;
        }
    }
}

// 将数据写入文件
void writeToFile(const vector<complex<double>> &Original_signal, const vector<complex<double>> &Proceeded_Signal, const string &filename, int Type_sel)
{
    ofstream file(filename);
    if (Type_sel==SEL_FFTIFFT)
    {
        file << "time(s)" << ", " << "Original Signal"<< ", " << "Filtered Signal" << "\n";
        for (size_t i = 0; i < Original_signal.size(); i++)
        {
            file << (double)i / SAMPLE_RATE << ", " <<Original_signal[i].real() << ", "<< Proceeded_Signal[i].real() << "\n";
       }
    }
    
    file.close();
}

// 测试模拟程序
int main()
{
    //FIR filter：50Hz选频FIR滤波器
    vector<double> 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<complex<double>> H(b.size());

    for(size_t i=0; i< b.size(); i++)
    {
        H[i] = complex<double>(b[i],0);
    }
    H.resize(N, complex<double>(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<complex<double>>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<vector<complex<double>>> 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<double>(0.0, 0.0));
        // devide Original_signal into 3 parts to simulate consequent DAQ for 3 cycles
        for(size_t j=0; j<DataSeg_len_L;j++)
        {
            seg_Dates[i][j] = Original_signal[i * DataSeg_len_L + j];
        }

        seg_Dates[i].resize(N, complex<double>(0.0, 0.0));                     // resize each data segment to N = L + M - 1
    }

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

    // (4-1) Simulate 3 cycles of involution: 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;
}

时间仓促，有些功能和细节并没有考虑太多，这里功能验证是第一。

FFT/IFFT是基于库里-图基的算法，各位可以选用其他的来优化替代；

有些参数可以单独变，有些参数却是关联的。

[1]《数字信号处理—原理、算法与应用》，第四版，普鲁克斯著，方艳梅等译，北京电子工业出版社，2014.8