【信号滤波 (下)】其他滤波器(一) 基于FFT的频域数字滤波

基于FFT的频域数字滤波原理

傅里叶变换和滤波算法去除ADC采样中的噪声
在之前介绍了对整个信号用傅里叶变换的方法看频谱，不用二阶滤波器的情况下，也可以直接在频谱上去除噪声所在的频率分量，再通过傅里叶逆变换(IFFT)的方法还原成时域信号。
该方法直接使用了傅里叶变换和傅里叶逆变换，并且直接修改频谱上噪声分量显得非常直观易懂，但是无法用于实时滤波。

滤波原理详解

一般步骤如下：

1. 使用快速傅里叶变换 （FFT） 将信号转换为频域；
2. 识别并设置噪声处陷波频率分量；
3. 将对应陷波频率的分量归零或减小；
4. 使用快速傅里叶逆变换（IFFT） 将信号转换回时域。

步骤详解

1. 实际应用中ADC采样不是连续的，而是每隔几毫米读一次ADC数据，获得是离散的时域信号，数学上用$x[n]$表示，n表示第几次采样。
   离散的时域$x[n]$转成离散频域$X[k]$的公式如下：
   $$X[k]=\sum _{n=0}^{N-1}x[n]\cdot e^{-i2\pi \frac{k}{N}n}$$
   其中:
   • $N$ 是采样总次数
   • $k$ 是频率的索引

上面这个公式人工算起来很麻烦，要一个个算过去，但计算机两个个for循环就解决了，大循环是k从0到N-1（因为再往后就循环了，没意义），小循环是n从0到N-1。

2. 傅里叶变换的频率分辨率取决于采样率和傅里叶变换索引的数量（信号的长度）。公式如下：
   $$f_k=\frac{k\cdot f_s}{N}$$
   其中:

   • $f_k$ 对应于第k个索引处频率
   • $f_s$ 是采样频率
   • $N$ 是采样总数 (FFT 长度).

3. 把频域中噪声处的频率归零，简单而且直观。

4. 傅里叶逆变换也是两个for循环就解决了，公式如下：
   $$x[n]=\frac{1}{N}\sum _{k=0}^{N-1}X[k]\cdot e^{i2\pi \frac{k}{N}n}$$

这个方法的缺陷是无法实时滤波，但对于非实时信号，比如说采集一个窗口时间内的信号后再分析，就可以用基于FFT的频域数字滤波。

matlab实现滤波算法

这里将基于FFT的频域数字滤波与滑动平均滤波和二阶陷波滤波一起做比较，写在matlab代码里，如下：

data = xlsread('importData.xlsx'); % 读取整个数据文件
voltage = data(1:500, 1); 

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

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

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);

ignore_below_frequency = 5;

% 滤除低频直流分量干扰
validIndices = f_used >= ignore_below_frequency; 
filteredP1 = P1(validIndices); 
filteredF = f_used(validIndices); 

% 寻找峰值
[peakAmplitude, peakIndex] = max(filteredP1); 
peakFrequency = filteredF(peakIndex); 

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

% 滑动平均滤波
window_size = fix(Fs/peakFrequency);
if mod(window_size, 2) == 0
    window_size = window_size + 1; % Ensure window size is odd for symmetry
end
filtered_voltage_MA = movmean(voltage, window_size);

subplot(3,2,3);
plot(time, filtered_voltage_MA, 'b');
xlabel('采样时间 (ms)');
ylabel('电压 (mV)');
title(['滑动平均滤波 窗口大小', num2str(window_size), ' 处理后']);
grid on;

% 设计二阶滤波器
wo = peakFrequency / (Fs / 2);          
bw = wo / 5;                 
notchFilter = designfilt('bandstopiir', ...
    'FilterOrder', 2, ...
    'HalfPowerFrequency1', wo - bw/2, ...
    'HalfPowerFrequency2', wo + bw/2, ...
    'DesignMethod', 'butter');

% 应用二阶滤波器
filtered_voltage = filtfilt(notchFilter, voltage); % Zero-phase filtering

subplot(3,2,4);
plot(time, filtered_voltage, 'm');
xlabel('采样时间 (ms)');
ylabel('电压 (mV)');
title(['陷波滤波 频率', num2str(peakFrequency), 'Hz 带宽', num2str(bw)]);
grid on;

k_noise = round(peakFrequency * N / Fs);  % 找到噪声频率

Y(k_noise) = 0;  % 实部归零
Y(k_noise + 1) = 0;  % 实部归零
Y(k_noise - 1) = 0;  % 实部归零
Y(N-k_noise + 1) = 0;  % 虚部归零
Y(N-k_noise) = 0;  % 虚部归零
Y(N-k_noise + 2) = 0;  % 虚部归零

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);

subplot(3,2,5);
plot(f_used, P1); % 绘制单边频谱的幅值谱
xlim([ignore_below_frequency 200]);
xlabel('频率 (Hz)'); % 频率轴标签（需要知道采样率才能正确标注）
ylabel('幅值');
title('去除噪声频点后的频域波形'); % 图表标题
grid on; % 显示网格

filtered_voltage_ifft = ifft(Y);  % 傅里叶逆变换
filtered_voltage_ifft_real = real(filtered_voltage_ifft);

subplot(3,2,6);
plot(time, filtered_voltage_ifft_real);
title('去除噪声频点后逆变换到时域的波形');
xlabel('采样时间 (ms)');
ylabel('电压 (mV)');
grid on;

基于FFT滤波效果对比

我这里对一段已经采集完成的时域信号滤波，在叠加20Hz、50Hz、80Hz频率噪声的情况下，基于FFT的频域数字滤波和滑动平均滤波以及陷波滤波相比，效果比后两者都好。



但缺点就是只能用在非实时滤波中，实时滤波除了前面提到的陷波滤波还有切比雪夫滤波器/卡尔曼滤波器等。

QT C++实现基于FFT的频域数字滤波

因为滤波一般用在嵌入式上，主要可用分成三个函数来写，傅里叶正变换函数，去噪归零函数，傅里叶逆变换函数：

QVector<double> FFT(const QVector<double>& voltageData) {
    int N = voltageData.size();
    
    // 创建复数形式的FFTW数组
    fftw_complex* in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    fftw_complex* out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    
    // 将输入数据转换为复数形式（仅限实部）
    for (int i = 0; i < N; ++i) {
        in[i][0] = voltageData[i];  // 实部
        in[i][1] = 0.0;             // 虚部
    }
    
    // 创建傅里叶正变换
    fftw_plan plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    
    // 执行
    fftw_execute(plan);
    
    // 输出包含复频域信号
    QVector<double> freqData(N);
    for (int i = 0; i < N; ++i) {
        freqData[i] = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]);  // Magnitude
    }

    fftw_destroy_plan(plan);
    fftw_free(in);
    fftw_free(out);
    
    return freqData;

去噪归零函数

void NotchToFrequencyZero(QVector<fftw_complex>& freqData, int notchFreq, double sampleRate) {
    int N = freqData.size();
    double freqResolution = sampleRate / N;
    
    // 寻找噪声频率
    int notchBin = static_cast<int>(notchFreq / freqResolution);
    
    // 噪声处频率分量归零
    int notchWidth = 1; // 设置归零带宽
    for (int i = notchBin - notchWidth; i <= notchBin + notchWidth; ++i) {
        if (i >= 0 && i < N) {
            freqData[i][0] = 0.0;  // 实部归零
            freqData[i][1] = 0.0;  // 虚部归零
        }
    }
}

傅里叶逆变换函数

QVector<double> IFFT(const QVector<fftw_complex>& freqData) {
    int N = freqData.size();
    
    fftw_complex* in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    fftw_complex* out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    

    for (int i = 0; i < N; ++i) {
        in[i][0] = freqData[i][0];  
        in[i][1] = freqData[i][1];  
    }
    
    // 创建傅里叶逆变换
    fftw_plan plan = fftw_plan_dft_1d(N, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
    
    // 执行
    fftw_execute(plan);
    
    // 输出包含时域信号
    QVector<double> timeDomainData(N);
    for (int i = 0; i < N; ++i) {
        timeDomainData[i] = out[i][0] / N;  // Normalize the real part
    }

    fftw_destroy_plan(plan);
    fftw_free(in);
    fftw_free(out);
    
    return timeDomainData;
}

基于FFT滤波的优势和劣势

首先说优势，通过和二阶陷波滤波对比，不会在起始和结束处抖动，原因也很好理解。我们看IIR滤波的时域下公式：
$$y[n]+a_{1}y[n-1]+a_{2}y[n-2]=b_{0}x[n]+b_{1}x[n-1]+b_{2}x[n-2]$$
n代表当前采样时间项，n-1 n-2代表前一两次采样时间项，那么在边缘端会缺失，导致振荡，而基于FFT滤波完全是基于频域和时域的转换，所以不会存在振荡。

当然劣势也很明显，不能实时滤波。这点也很好理解，我们知道傅里叶变换最早的要求是信号是周期重复的才可以，离散傅里叶变换把不是周期的信号，通过重复的方法延申成周期信号，所以必须要又一个时间段来采集信号才能做傅里叶变换。
而陷波滤波器的实现方式可以电路分析里的频率响应对应上，可以理解成RC滤波电路/LC滤波电路的数字化，而RC/LC滤波电路在可以把特定频率的信号衰减，不需要采样时间窗口，所以可以实现实时滤波。