FFT变化与反变换滤波分析

using System;
using MathNet.Numerics;
using MathNet.Numerics.IntegralTransforms;
using MathNet.Numerics.LinearAlgebra;

public class Filter
{
    public static double[] Apply(double[] input, int dimension, double lowerCutoff, double higherCutoff, double sampleRate)
    {
        // 计算 FFT
        Complex32[] fftInput = Array.ConvertAll(input, x => new Complex32(Convert.ToSingle(x), 0));
        Fourier.Forward(fftInput, FourierOptions.Matlab);

        // 创建掩码
        int n = fftInput.Length;
        double[] frequencies = new double[n];
        for (int i = 0; i < n; i++)
        {
            frequencies[i] = i * (sampleRate / n);
        }
        bool[] mask = new bool[n];
        for (int i = 0; i < n; i++)
        {
            mask[i] = frequencies[i] > lowerCutoff && frequencies[i] < higherCutoff;
        }

        // 应用掩码
        for (int i = 0; i < n; i++)
        {
            if (!mask[i])
            {
                fftInput[i] = Complex32.Zero;
            }
        }

        // 计算 IFFT
        Fourier.Inverse(fftInput, FourierOptions.Matlab);

        // 转换回 double 数组
        double[] output = Array.ConvertAll(fftInput, x => Convert.ToDouble(x.Real));

        return output;
    }
}

以下为对应Matlab代码，转载自他人博客

% FILTERED = ideal_bandpassing(INPUT,DIM,WL,WH,SAMPLINGRATE)
% 
% Apply ideal band pass filter on INPUT along dimension DIM.
% 
% WL: lower cutoff frequency of ideal band pass filter
% WH: higher cutoff frequency of ideal band pass filter
% SAMPLINGRATE: sampling rate of INPUT
% 
% Copyright (c) 2011-2012 Massachusetts Institute of Technology, 
% Quanta Research Cambridge, Inc.
%
% Authors: Hao-yu Wu, Michael Rubinstein, Eugene Shih, 
% License: Please refer to the LICENCE file
% Date: June 2012
%
function filtered = ideal_bandpassing(input, dim, wl, wh, samplingRate)
 
    if (dim > size(size(input),2))
        error('Exceed maximum dimension');
    end
 
    input_shifted = shiftdim(input,dim-1);
    Dimensions = size(input_shifted);
    
    n = Dimensions(1);
    dn = size(Dimensions,2);
    
    
    Freq = 1:n;
    Freq = (Freq-1)/n*samplingRate;
    mask = Freq > wl & Freq < wh;
    
    Dimensions(1) = 1;
    mask = mask(:);
    mask = repmat(mask, Dimensions);
 
    
    F = fft(input_shifted,[],1);
    
    F(~mask) = 0;
    
    filtered = real(ifft(F,[],1));
    
    filtered = shiftdim(filtered,dn-(dim-1));
    
end