【信号滤波 (上)】傅里叶变换和滤波算法去除ADC采样中的噪声(Matlab/C++)

一、ADC采样的噪声简介

经常用ADC采样的朋友都知道，采样的数据不是纯净的信号，会伴随着干扰。
一般把周期性的干扰叫噪声，既然是周期性的，一定能算出周期时长T，取倒数就是频率f，单位Hz。

1.1 常见的ADC噪声来源

电源是噪声的一大来源，如果用的是交流电，交流电通过整流电路中的二极管元件，过滤掉交流电的负半周，然后通过滤波电容滤波让交流电趋于平缓，形成直流电。交流信号肯定不能滤干净的，这就是电源纹波产生的一大原因，电源纹波会导致ADC采样信号出现噪声。

二、信号的时域到频域转换

刚接触信号的时候搞不清楚什么是时域什么是频域，大学里的《信号与系统》《控制工程与原理》《信号处理》等书上提过，但因为没有结合实际使用学过都忘记了。
时域可以理解成横轴是时间t，纵轴是其他有幅值变换的物理量(如电压，声响，亮度，温度，振幅等)，如下图左图所示；
频域是横轴为频率大小，纵轴为频率分量的强度(或者叫功率密度)；
而傅里叶变换就是时域和频域的联结。

2.1 傅里叶变换

从时域变换到频域用的是傅里叶变换，傅里叶觉得计算波形相加和分解太麻烦了，发明了一种简单的方法让波形相加能像“1+1=2”一样直观，他的方法就是把波形从时域转换到频域再相加或者分解，计算好后再转回时域波形。我们把傅里叶发明的这套方法叫做“傅里叶变换”和“傅里叶逆变换”。
公式如下，如果看不懂一是回去复习微积分，二是用我下面的方法巧记：
F ( ω ) = F { f ( t ) } = 1 / 2 π ∫ − ∞ ∞ f ( t ) e − i ω t   d t ( 傅里叶变换 ) F(\omega) = \mathcal{F}\{f(t)\} ={1/2\pi} \int_{-\infty}^{\infty} f(t) e^{-i \omega t} \, dt (傅里叶变换) F(ω)=F{f(t)}=1/2π∫−∞∞​f(t)e−iωtdt(傅里叶变换)
f ( t ) = F − 1 { F ( ω ) } = ∫ − ∞ ∞ F ( ω ) e i ω t   d ω ( 傅里叶逆变换 ) f(t) = \mathcal{F}^{-1}\{F(\omega)\} = \int_{-\infty}^{\infty} F(\omega) e^{i \omega t} \, d\omega (傅里叶逆变换) f(t)=F−1{F(ω)}=∫−∞∞​F(ω)eiωtdω(傅里叶逆变换)
也可以用$\omega =2\pi f$带入
F ( f ) = F { f ( t ) } = ∫ − ∞ ∞ f ( t ) e − 2 π i f t   d t F(f) = \mathcal{F}\{f(t)\} = \int_{-\infty}^{\infty} f(t) e^{-2\pi i f t} \, dt F(f)=F{f(t)}=∫−∞∞​f(t)e−2πiftdt
f ( t ) = F − 1 { F ( f ) } = ∫ − ∞ ∞ F ( f ) e 2 π i f t   d f f(t) = \mathcal{F}^{-1}\{F(f)\} = \int_{-\infty}^{\infty} F(f) e^{2\pi i f t} \, df f(t)=F−1{F(f)}=∫−∞∞​F(f)e2πiftdf

巧记傅里叶变换

先看这个动画
傅里叶变换动画
首先看懂公式最里面的这个$e^{-i\omega t}$，这其实就是三角函数在复数域的延申。
已知欧拉公式$e^{-i\omega t}=cos(\omega t)-isin(\omega t)$，$e^{-i\omega t}$可以看成一个距离原地1单位的球在绕着原点以$\omega$角速度旋转，从x轴(实数)看过去就是$cos(\omega t)$，从y轴(虚数)看过去就是$sin(\omega t)$。

然后再看$f(t)e^{-i\omega t}$，下面动图最上面的蓝线是时域函数f(t)，红线就是$e^{-i\omega t}$的实数轴部分，这里有两个变量$\omega$和t把大家看懵了，我们细细捋一下，下面的横轴就是t，t从(0,∞)的变换都已经显示在数轴上了，可以看成静态量，导致红线变化的是另一个变量$\omega$，它才是动态变量。

$f(t)e^{-i\omega t}$也就是把$f(t)$和$e^{-i\omega t}$相乘，也就是中间图。可以看到当$e^{-i\omega t}$变化的和$f(t)$波动频率一致的时候，乘积结果就是$si{n}^2(t)=(1-cos(2t))/2$。

最后经过积分，${\int }_{-\infty }^{\infty }f(t)e^{-i\omega t}dt$得到的下图。注意看，这是对时间t进行积分，之前说过时间是静态量，所以图中横轴已经不是时间t了，而是角频率$\omega$，公式$F(\omega )=1/2\pi {\int }_{-\infty }^{\infty }f(t)e^{-i\omega t}dt$表示的是积分结果随着频率变化而变化。

最后一张图和时间已经没关系了，横轴已经是频率了，这样就完成了时域到频域的转化。

PS:为啥要费劲心思把周期信号从时域转频域，关键就是"周期"，周期意味着经过一定时间后会重复，所以时间反而不重要了，重要的是频率。下图的最高点对应的就是原始信号的频率。

(动图参考来源：https://blog.csdn.net/YuHWEI/article/details/136678900)

上面的图能解释傅里叶变换但不是真正的傅里叶变换，因为我们只计算了实数轴，真正的傅里叶变换后，正弦信号是在8.3Hz上的一道冲击信号。

同理，多个周期性正弦波叠加也可以从频域上看出频率分量。
在频域上做周期性信号的合成与分解是不是就简单多了呢。

通过不同频率的正弦波叠加，想要什么波形都能叠加出来。
反过来，任何一种波形，都可以分解成不同频率的正弦波叠加。

三、傅里叶变换和滤波算法工程实现

3.1 使用Matlab计算信号时域到频域的变换

matlab使用快速傅里叶变换函数(fft)，注意fft后的复数域的量，还需要abs()一下才能显示到实数轴上。

data = xlsread('origin_signal_data.xlsx'); % 读取原始波形数据
voltage = data(:, 1); 

N = length(voltage); % 数据长度

t = 2.5; %总采样时间
T = t/N; %采样周期=总采样时间/总采样次数
Fs = 1/T; %采样率是1s内的采样次数
time = 0:5:2500-5; %设置横坐标数组
plot(time, voltage, 'r', 'LineWidth', 2); % 绘制单边频谱的幅值谱
xlabel('采样时间 (ms)'); % 频率轴标签（需要知道采样率才能正确标注）
ylabel('电压 (mV)'); % 幅值轴标签
title('原始时域波形'); % 图表标题
grid on;

% data_no_dc = timeSeries - signal_mean;
Y = fft(voltage); % 计算FFT
P2 = abs(Y/N); % 双边频谱的幅值，并归一化
P1 = P2(1:N/2+1); % 单边频谱的幅值（正频率部分）
P1(2:end-1) = 2*P1(2:end-1); % 因为FFT结果是对称的，所以正频率部分的幅值要乘以2（除了第一个和最后一个点）

f = (0:N-1)*(1/N); % 频率向量（归一化频率）
f_actual = f*(Fs); % 实际频率（需要知道采样率）
f_used = f_actual(1:N/2+1); %使用到的只有一半

plot(f_used, P1); % 绘制单边频谱的幅值谱
xlabel('频率 (Hz)'); % 频率轴标签（需要知道采样率才能正确标注）
ylabel('幅值');
title('傅里叶变换后频域波形'); % 图表标题
grid on; % 显示网格

下面左图的原始数据中还有周期性噪声干扰，通过傅里叶变换后看到噪声主要集中在60Hz。

注意！我这里是去除了0Hz的数据，0Hz是信号中的直流分量，信号都是直流分量和交流分量的叠加，如果傅里叶变换后0Hz有很高的幅值，那就是直流分量，我们要处理的是噪声肯定是交流分量，所以作图时不要把0Hz算进去。

3.2 使用Matlab去除特定频点噪声

寻找峰值算噪声频率

ignore_below_frequency = 1; % 1Hz以下不管，因为是信号的直流分量

% 寻找傅里叶变换后最大幅值所在的频点
validIndices = f_used >= ignore_below_frequency; 
filteredP1 = P1(validIndices); 
filteredF = f_used(validIndices); 

[peakAmplitude, peakIndex] = max(filteredP1); 
peakFrequency = filteredF(peakIndex); % 获得噪声频率

构建陷波滤波器滤除噪声频点

构建陷波滤波器，设置滤波频点和带宽

% 设置60Hz陷波滤波器
wo = peakFrequency / (Fs / 2);           
bw = wo / 20;                 % 设置滤波器带宽
notchFilter = designfilt('bandstopiir', ...
    'FilterOrder', 2, ...
    'HalfPowerFrequency1', wo - bw/2, ...
    'HalfPowerFrequency2', wo + bw/2, ...
    'DesignMethod', 'butter');

% 使用滤波器
filtered_voltage = filtfilt(notchFilter, voltage); 

陷波滤波器与滑动平均滤波器对比

滑动平均滤波器就是设置一个窗口，窗口内包含多个采样频点，比如前7个数据为一个窗口，取窗口内数据的平均值为当前点的值，然后窗口向后滑动，滑动平均算法的核心思想是在一个正弦周期内所有点之和为零。
如果选的窗口刚好是噪声周期，那就滤除干扰了。
但滑动平均在较低频率下好用，高频后窗口大小无法继续缩小，不如陷波滤波。

3.3 C++ 实现ADC采样信号滤波

matlab验证方法，工程实现还是要用C++，这里我要把算法最终应用到嵌入式设备上去，所以使用QT的C++来部署验证。

QT C++需要使用快速傅里叶变换库和矩阵库支持
在.pro文件里加入fftw3库。

LIBS += -lfftw3

同时下载Eigen库放到项目文件中，并加入头文件路径

INCLUDEPATH += ./eigen-3.4.0

陷波滤波器没有现场的库，需要手动写

QVector<double> applyNotchFilter(const QVector<double> &voltageData, double sampleRate, double notchFreq, double qualityFactor)
{
    QVector<double> filteredData(voltageData.size());
    double w0 = 2.0 * M_PI * notchFreq / sampleRate;
    double alpha = sin(w0) / (2.0 * qualityFactor);
    double b0 = 1.0;
    double b1 = -2.0 * cos(w0);
    double b2 = 1.0;
    double a0 = 1.0 + alpha;
    double a1 = -2.0 * cos(w0);
    double a2 = 1.0 - alpha;
    b0 /= a0;
    b1 /= a0;
    b2 /= a0;
    a1 /= a0;
    a2 /= a0;
    double x1 = 0.0, x2 = 0.0; 
    double y1 = 0.0, y2 = 0.0; 
    for (int i = 0; i < voltageData.size(); ++i) {
        double x0 = voltageData[i];
        double y0 = b0 * x0 + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2;
        filteredData[i] = y0;
        x2 = x1;
        x1 = x0;
        y2 = y1;
        y1 = y0;
    }
    return filteredData;
}

为了像matlab一样能直观的看到滤波效果，用QChart作图

void displayFilterResult(const QVector<double> timeData, const QVector<double> &voltageData, const QVector<double> filteredData, const double maxFrequency)
{
    // 创建QLineSeries对象，并设置它们的名称
    QLineSeries *seriesVoltage = new QLineSeries();
    seriesVoltage->setName("原始电压数据"); // 设置名称为“原始电压数据”
    QLineSeries *seriesFiltered = new QLineSeries();
    seriesFiltered->setName("滤波后电压数据"); // 设置名称为“滤波后电压数据”

    // 填充数据到QLineSeries对象
    for (int i = 0; i < timeData.size(); ++i) {
        seriesVoltage->append(timeData[i], voltageData[i]);
        seriesFiltered->append(timeData[i], filteredData[i]);
    }

    // 创建QChart对象，并添加系列（与之前相同）
    QChart *chart = new QChart();
    chart->addSeries(seriesVoltage);
    chart->addSeries(seriesFiltered);

    // 默认情况下，图例是可见的，但你可以通过以下方式访问和配置它
    QLegend *legend = chart->legend();
    legend->setVisible(true); // 通常情况下这是默认的，但你可以显式设置
    legend->setAlignment(Qt::AlignBottom); // 设置图例的位置，这里是在底部

    // 为图表设置标题（与之前相同）
    QString frequency = QString::number(maxFrequency, 'f', 1) + "Hz";
    chart->setTitle("滤波前后对比——干扰频率" + frequency);

    // 创建QValueAxis对象并设置轴的范围
    QValueAxis *axisX = new QValueAxis();
    axisX->setLabelFormat("%g s");
    chart->addAxis(axisX, Qt::AlignBottom);
    seriesVoltage->attachAxis(axisX);
    seriesFiltered->attachAxis(axisX);

    QValueAxis *axisY = new QValueAxis();
    axisY->setLabelFormat("%g V");
    chart->addAxis(axisY, Qt::AlignLeft);
    seriesVoltage->attachAxis(axisY);
    seriesFiltered->attachAxis(axisY);

    // 为不同的系列设置不同的颜色
    seriesVoltage->setPen(QPen(QColor(255, 0, 0))); // 红色
    seriesFiltered->setPen(QPen(QColor(0, 0, 255))); // 蓝色

    // 创建QChartView对象并设置QChart
    QChartView *chartView = new QChartView(chart);
    chartView->setRenderHint(QPainter::Antialiasing);

    // 创建主窗口并显示图表
    QMainWindow* window = new QMainWindow();
    window->setCentralWidget(chartView);
    window->resize(1920, 1080);
    window->show();
}

主函数如下

#include <QFile>
#include <QString>
#include <QVector>
#include <QDebug>
#include <QtCharts/QChartView>
#include <QtCharts/QLineSeries>
#include <QtCharts/QValueAxis>
#include <QMainWindow>
#include <QApplication>

#include <fftw3.h>
#include <Eigen/Dense>


using namespace Eigen;
using namespace std;

QT_CHARTS_USE_NAMESPACE

int main(int argc, char *argv[]) {
    QApplication a(argc, argv);

    QFile inputFile("./source_signal_data.csv");
    if (!inputFile.open(QIODevice::ReadOnly | QIODevice::Text)) {
        qDebug() << "文件打开失败";
        return -1;
    }

    QVector<double> voltageData;
    while (!inputFile.atEnd()) {
        QByteArray line = inputFile.readLine();
        voltageData.push_back(line.toDouble());
    }
    inputFile.close();

    // 配置周期和频率
    int N = voltageData.size();
    double t = 2.5; // 总时间
    double T = t / N;
    double Fs = 1 / T; // 采样频率

    // 生成时间数据，从0s开始，每隔0.005s到2.5s
    QVector<double> timeData;
    for (double t = 0; t < (2.5 - 0.005); t += 0.005) {
        timeData.append(t);
    }

    // 将信号数据放到实部中
    VectorXd voltageEigen(N);
    for (int i = 0; i < N; ++i) {
        voltageEigen[i] = voltageData[i];
    }

    fftw_complex* in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
    fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
    fftw_plan plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

    for (int i = 0; i < N; i++) {
        in[i][0] = voltageEigen[i]; // 实部
        in[i][1] = 0.0;             // 虚部
    }

    fftw_execute(plan);  //应用傅里叶变换

    VectorXcd fftResult(N);
    for (int i = 0; i < N; i++) {
        fftResult[i] = std::complex<double>(out[i][0], out[i][1]);
    }

    // 清楚缓存
    fftw_destroy_plan(plan);
    fftw_free(in);
    fftw_free(out);

    // 制作频谱
    VectorXd absFFT = fftResult.cwiseAbs() / N;
    VectorXd P1 = absFFT.head(N / 2 + 1);
    P1.segment(1, P1.size() - 2) *= 2;

    VectorXd f_actual = VectorXd::LinSpaced(N / 2 + 1, 0, Fs / 2);

    // 滤除直流分量
    QVector<double> filteredP1;
    QVector<double> filteredF;
    for (int i = 0; i < f_actual.size(); ++i) {
        if (f_actual(i) >= 1) {
            filteredP1.push_back(P1(i));
            filteredF.push_back(f_actual(i));
        }
    }

    // 导出到csv
    QFile outputFile("../FFT_test/filtered_data.csv");
    if (outputFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
        QTextStream out(&outputFile);
        out << "频率,幅值\n";
        for (int i = 0; i < filteredF.size(); ++i) {
            out << filteredF[i] << "," << filteredP1[i] << "\n";
        }
        outputFile.close();
    } else {
        qDebug() << "导出失败" ;
    }

    // 寻找噪声最大的幅值和对应的频率
    auto maxIt = max_element(filteredP1.begin(), filteredP1.end());
    double maxAmplitude = 0.0;
    double maxFrequency = 0.0;
    if (maxIt != filteredP1.end())
    {
        size_t maxIndex = distance(filteredP1.begin(), maxIt);
        maxAmplitude = *maxIt;
        maxFrequency = filteredF[maxIndex];
    }
    else
    {
        qDebug() << "未找到噪声频率";
    }

    double qualityFactor = 1; // 设置陷波滤波带宽
    QVector<double> filteredData = applyNotchFilter(voltageData, Fs, maxFrequency, qualityFactor);

    // 导出到csv
    QString notchFileName = "./notch_filtered_" + QString::number(maxFrequency, 'f', 1) + "Hz.csv";
    QFile notchOutputFile(notchFileName);
    if (notchOutputFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
        QTextStream out(&notchOutputFile);
        for (int i = 0; i < filteredData.size(); ++i) {
            out << filteredData[i] << "\n";
        }
        notchOutputFile.close();
    } else {
        qDebug() << "导出失败";
    }

    //在图表上显示滤波前后效果
    displayFilterResult(timeData, voltageData, filteredData, maxFrequency);

    return a.exec();
}

QT滤波实现效果如下：