java通过频域进行低通滤波

滤波分为2种，一种是时域域滤（比如Iir数字滤波器），一种是频域滤波（先傅里叶正变换生成频域数据，在频域将不需要的数据置零，然后反傅里叶变换生成时域曲线）。本文采用的是频域低通滤波。

 

/**
     * 频域滤波
     * @param originalPoints 原始数据
     * @param samplingInterval 采样间隔, 单位是秒(s)
     * @param lowpass 低通
     * @return 滤波之后的数据
     */
    public static List<SamplingPoint> lowpassFromFrequencyDomain(List<SamplingPoint> originalPoints, double samplingInterval, int lowpass){
        int dataLen = originalPoints.size();
        double[] signal = new double[dataLen];
        for(int i = 0; i < originalPoints.size(); i++){
            signal[i] = originalPoints.get(i).getRealY();
        }

        double Fs = 1.0 / samplingInterval;//采样频率
        double frequencyResolution = Fs / dataLen;//频率分辨率

        //进行傅里叶变换,使用的是commons-math包里面的快速傅里叶算法
        FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);
        Complex[] complexArr = fft.transform(signal, TransformType.FORWARD);//正向傅里叶变换

        //进行滤波,注意条件
        if(lowpass < Fs / 2.0){
            for(int i = 0; i < complexArr.length; i++){
                //对称
                if(i * frequencyResolution > lowpass && i * frequencyResolution < (Fs - lowpass)){
                    complexArr[i] = new Complex(0, 0);
                }
            }
        }

        //反傅里叶变换
        Complex[] timeDomainArr = fft.transform(complexArr, TransformType.INVERSE);
        List<SamplingPoint> points = new ArrayList<>();
        for(int i = 0; i < timeDomainArr.length; i++){
            //只用获取实部，不用获取虚部(虚部理论上应该是0),实部的数据就是时域曲线
            SamplingPoint point = new SamplingPoint(i, (float) timeDomainArr[i].getReal());
            points.add(point);
        }
        return points;
    }

代码说明：SamplingPoint是一个javabean，封装了实测数据和一些其他数据，getRealY()获取到的就是实测数据。这个函数中使用了commons-math.jar包，commons-math是apache提供的一个数学计算包。如果不希望引入这个jar包，可以使用下面这个函数进行傅里叶正变换和反变换，下面的fft()函数是c语言代码，可以自己转换成java代码。

/**********************************************************************
  FFT - calculates the discrete fourier transform of an array of double
  precision complex numbers using the FFT algorithm.

  c = pointer to an array of size 2*N that contains the real and
    imaginary parts of the complex numbers. The even numbered indices contain
    the real parts and the odd numbered indices contain the imaginary parts.
      c[2*k] = real part of kth data point.
      c[2*k+1] = imaginary part of kth data point.
  N = number of data points. The array, c, should contain 2*N elements
  isign = 1 for forward FFT, -1 for inverse FFT.
*/

void FFT( double *c, int N, int isign )
{
  int n, n2, nb, j, k, i0, i1;
  double wr, wi, wrk, wik;
  double d, dr, di, d0r, d0i, d1r, d1i;
  double *cp;

  j = 0;
  n2 = N / 2;
  for( k = 0; k < N; ++k )
  {
    if( k < j )
    {
      i0 = k << 1;
      i1 = j << 1;
      dr = c[i0];
      di = c[i0+1];
      c[i0] = c[i1];
      c[i0+1] = c[i1+1];
      c[i1] = dr;
      c[i1+1] = di;
    }
    n = N >> 1;
    while( (n >= 2) && (j >= n) )
    {
      j -= n;
	  n = n >> 1;
    }
    j += n;
  }

  for( n = 2; n <= N; n = n << 1 )
  {
    wr = cos( 2.0 * M_PI / n );
    wi = sin( 2.0 * M_PI / n );
    if( isign == 1 ) wi = -wi;
    cp = c;
    nb = N / n;
    n2 = n >> 1;
    for( j = 0; j < nb; ++j )
    {
      wrk = 1.0;
      wik = 0.0;
      for( k = 0; k < n2; ++k )
      {
        i0 = k << 1;
        i1 = i0 + n;
        d0r = cp[i0];
        d0i = cp[i0+1];
        d1r = cp[i1];
        d1i = cp[i1+1];
        dr = wrk * d1r - wik * d1i;
        di = wrk * d1i + wik * d1r;
        cp[i0] = d0r + dr;
        cp[i0+1] = d0i + di;
        cp[i1] = d0r - dr;
        cp[i1+1] = d0i - di;
        d = wrk;
        wrk = wr * wrk - wi * wik;
        wik = wr * wik + wi * d;
      }
      cp += n << 1;
    }
  }
}