FFTフィルタ各種 – C言語でFFTWでMATLABで

過去に自分が作ったFFTフィルタをいくつか見つけた.

動作が保証出来るのはMATLABのとC言語のかな.
FFTWを使ったやつは・・なんでか無理にC++で書いてるもんだから,
微妙やもしれない(結局使わずなので,動作確認して無い).
まあ,多分動くんじゃないかなーとは思うけども(何だそれ).

画像処理のための(2次元)FFTフィルタ~~みたいな記事はよく見かけるのに,
何故か(生体)信号処理やらの話ではFIRフィルタばっかりで,
FFTフィルタの話が書かれていないんだよね.
(厳密には,「FFTフィルタ」じゃなくて「FFTを利用したデジタル・フィルタ」か)

FIRフィルタに比べてFFTフィルタを用いたほうが,
位相遅れが無いので元信号に対して位相特性を考慮しなくて良いし,
FIRフィルタに比べて遥かに急峻(欲しい部分以外0に落とすから)なフィルタが得られるから,
信号処理においてはFIRフィルタよりもFFTフィルタのほうがメリット多そうなんだけど.

まあ,それはさておき・・・・・・
以下,MATLAB,FFTW,C言語で書いたFFTフィルタ.
ちなみに,どれも計算結果はノーマライズしてる.

————————————————————–

            MATLABでフーリエ・デジタル・フィルタ

————————————————————–
FFTフィルタ(いちおう多次元対応のつもり;1次元しか確認して無いが)

%FFT filter
%fftfilter(Sn, sf, lf, hf)
% Sn:  時系列データ
% sf: サンプリングレート
% lf:  ローカット周波数
% hf:  ハイカット周波数
%
function y=fftfilter(Sn, sf, lf, hf)

[Rows, Columns]=size(Sn);      % データ数の取得
d=round(Rows/sf);
if mod(d,2)==1, d=d+1; end     % データ数が奇数の場合は偶数に変更

lp=round(lf*d);                % 低域周波数に相当するのデータ地点
hp=round(hf*d);                % 高域周波数に相当するのデータ地点
F=zeros(Rows, Columns);        % フィルタ用配列を0で初期化
F(lp:hp,Columns)=1;            % 実数領域のフィルタ形状の定義
F(Rows-hp:Rows-lp,Columns)=1;  % 虚数領域のフィルタ形状の定義

S=fft(Sn);                     % ノイズの時間情報を周波数情報に変換
S=S.*F;                        % フィルタリング
S=ifft(S);                     % 周波数情報を時間情報に変換
S=real(S);
S=S/max(abs(S));               % -1から1で標準化
y=S;                           %fftfilter結果を返す

end

————————————————————–
            MATLABでフーリエ・デジタル・フィルタ その2
————————————————————–
FFTフィルタ+フィルタ後の振幅およびパワースペクトル
(自分用;急ぎで作ったものなので汚い)
(音声解析のためサンプリングはsf=44100で固定)

function y=NrsCllSpmRsl(Sn)

sf=44100;                  %サンプリング周波数
lf=2000;                    %ローカットオフ周波数
hf=5000;          %ハイカットオフ周波数

[Rows, Columns]=size(Sn);      % データ数の取得
d=round(Rows/sf);
if mod(d,2)==1, d=d+1; end     % データ数が奇数の場合は偶数に変更

lp=round(lf*d);                % 低域周波数に相当するのデータ地点
hp=round(hf*d);                % 高域周波数に相当するのデータ地点
F=zeros(Rows, Columns);        % フィルタ用配列を0で初期化
F(lp:hp,Columns)=1;            % 実数領域のフィルタ形状の定義
F(Rows-hp:Rows-lp,Columns)=1;  % 虚数領域のフィルタ形状の定義

S=fft(Sn);                     % ノイズの時間情報を周波数情報に変換
S=S.*F;                        % フィルタリング
S=ifft(S);                     % 周波数情報を時間情報に変換
S=real(S);
S=S/max(abs(S));               % -1から1で標準化
T2=S;                           %fftfilter結果を返す

y=T2;

t=1:Rows;
t=t/sf;
S=y;

format long         %倍精度で表示のため

n=length(S);

T=fft(S);
T=T.*conj(T);
T=T+0.000000001;
T=T./max(T(:));

y2=sqrt(T);
T=20*log10(T);

subplot(3,1,1);
plot(t,S);
%axis([1 d -1 1]);
xlabel(‘Time [s]’);
ylabel(‘Amplitude’);
set(gca,’FontName’,’MS P明朝’,’FontSize’,12);
h=get(gca,’title’);
set(h,’fontsize’,12);
h=get(gca,’xlabel’);
set(h,’fontsize’,12);
h=get(gca,’ylabel’);
set(h,’fontsize’,12);

subplot(3,1,2);
x=linspace(0,sf/2,n/2);
plot(x,y2(1:n/2));
axis([1 sf/2 0 1.1]);
xlabel(‘Frequency [Hz]’);
ylabel(‘Amplitude’);
set(gca,’FontName’,’MS P明朝’,’FontSize’,12);
h=get(gca,’title’);
set(h,’fontsize’,12);
h=get(gca,’xlabel’);
set(h,’fontsize’,12);
h=get(gca,’ylabel’);
set(h,’fontsize’,12);

subplot(3,1,3);
x=linspace(0,sf/2,n/2);
plot(x,T(1:n/2));
axis([1 sf/2 -150 20]);
xlabel(‘Frequency [Hz]’);
ylabel(‘Power [dB]’);
set(gca,’FontName’,’MS P明朝’,’FontSize’,12);
h=get(gca,’title’);
set(h,’fontsize’,12);
h=get(gca,’xlabel’);
set(h,’fontsize’,12);
h=get(gca,’ylabel’);
set(h,’fontsize’,12);

%format              %デフォルトに戻す

end

————————————————————–
          FFTWを使ってフーリエ・デジタル・フィルタ
————————————————————–
FFTWを使ったFFTフィルタ(動作未確認)

//フーリエ・デジタル・フィルター
//解析結果はそのまま**acclに返ってくる
void Fourier_Filter(int size, double **accl)
{
    double HPF;
    double LPF;
    fftw_complex *in;
    fftw_complex *out;
    fftw_plan p;

    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * size);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * size);
   
    if( !in || !out )
    {
        fprintf( stderr, "failed to allocate %d[byte] memory(-.-)n", size);
        return false;
    }
   
    //HPF, LPF(ロー,ハイカットオフ周波数)を読み込み
    cout << "HPF:"        cin >> HPF;
    cout << "LPF:"        cin >> LPF;

    //フーリエ変換を実行
    for(int j=0; j<NUM_AXIS ; j++)
    {
        p = fftw_plan_dft_1d( size, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
        for(int i=0; i<size; i++)
        {
            in[i][0] = accl[i][j]);
            in[i][1] = 0.0;
        }
        fftw_execute(p);
        for(int i=0; i<HPF ; i++)            out[i][0] = 0;
        for(int i=0; i<size-LPF ; i++)        out[i][1] = 0;
        for(int i=0; i>LPF ; i++)            out[i][0] = 0;
        for(int i=0; i>size-HPF ; i++)        out[i][1] = 0;

        for(int i=0; i<size ; i++)
        {
            in[i][0] = out[i][0];
            in[i][1] = out[i][1];
        }

        //逆フーリエ変換を実行
        p = fftw_plan_dft_1d( size, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
        fftw_execute(p);

        for(int i=0 ; i<size ; i++)
            accl[i][j] = out[i][0];
    }
   
    //すべて解放
    if( p   )    fftw_destroy_plan(p);
    if( in  )    fftw_free(in);
    if( out )    fftw_free(out);

    cout << "フーリエ・デジタル・フィルターによる処理完了.." << endl;
}

————————————————————–

          C言語でフーリエ・デジタル・フィルタ

————————————————————–

/*— FFTFilter.h —*/
#if _MSC_VER > 1000        //VC++4.0以降
#pragma once
#endif // _MSC_VER > 1000

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int FFTFilter(double, double, double)

static void FFT(double*, double*, int);
static void IFFT(double*, double*, int);
static int log2(int);
static int pow2(int);
static double max(const double*, int);
/*— end of FFTFilter.h —*/

/*— FFTFilter.cpp —*/

#include "FFTFilter.h"

void FFTFilter(double *TMData, double LFF, double HFF)
{
    int i;
    int N;
    int LffVl = round(LFF);        //ローカットオフ周波数;整数値に丸める
    int HffVl = round(HFF);        //ハイカットオフ周波数;整数値に丸める
    double *FFTRsl_rl;
    double *FFTRsl_img;

    N = sizeof(TMData)/sizeof(double);
    FFTRsl_rl  = calloc(N, sizeof(double)); //メモリの確保
    FFTRsl_img = calloc(N, sizeof(double)); //メモリの確保

    for(i=0 ; i<N ; i++)    FFTRsl_rl[i]  = TMData[i];

    FFT(FFTRsl_rl, FFTRsl_img, N); //FFTの計算結果はそれぞれ引数に

    for(i=0; i < LffVl     ; i++)        FFTRsl_rl[i]  = 0.0;
    for(i=0; i < N-HffVl ; i++)        FFTRsl_img[i] = 0.0;
    for(i=0; i > HffVl   ; i++)        FFTRsl_rl[i]  = 0.0;
    for(i=0; i > N-LffVl ; i++)        FFTRsl_img[i] = 0.0;

    IFFT(FFTRsl_rl, FFTRsl_img, N);    //IFFTの計算結果はそれぞれ引数に

    MaxVl = max(FFTRsl_rl, N);
    for(i=0 ; i<N ; i++)    TMData[i] = FFTRsl_rl[i]/MaxVl;

    free(FFTRsl_rl);    /* メモリの解放 */
    free(FFTRsl_img);    /* メモリの解放 */
}

static void FFT(double x_real[], double x_imag[], int N)
{
  int i;
  int j;
  int k;
  int n;
  int m;
  int r;
  int stage;
  int NumOfStage;
  int *Indx;
  double a_real;
  double a_imag;
  double b_real;
  double b_imag;
  double c_real;
  double c_imag;
  double real;
  double imag;
 
  /* FFTの段数 */
  NumOfStage = log2(N);
 
  /* バタフライ計算 */
  for (stage = 1; stage <= NumOfStage; stage++)
    for (i = 0; i < pow2(stage – 1); i++)
      for (j = 0; j < pow2(NumOfStage – stage); j++)
      {
        n = pow2(NumOfStage – stage + 1) * i + j;
        m = pow2(NumOfStage – stage) + n;
        r = pow2(stage – 1) * j;
        a_real = x_real[n];
        a_imag = x_imag[n];
        b_real = x_real[m];
        b_imag = x_imag[m];
        c_real = cos((2.0 * M_PI * r) / N);
        c_imag = -sin((2.0 * M_PI * r) / N);
        if (stage < NumOfStage)
        {
          x_real[n] = a_real + b_real;
          x_imag[n] = a_imag + b_imag;
          x_real[m] = (a_real – b_real) * c_real – (a_imag – b_imag) * c_imag;
          x_imag[m] = (a_imag – b_imag) * c_real + (a_real – b_real) * c_imag;
        }
        else
        {
          x_real[n] = a_real + b_real;
          x_imag[n] = a_imag + b_imag;
          x_real[m] = a_real – b_real;
          x_imag[m] = a_imag – b_imag;
        }
      }
 
  /* インデックスの並び替えのためのテーブルの作成 */
  Indx = calloc(N, sizeof(int));
  for (stage = 1; stage <= NumOfStage; stage++)
    for (i = 0; i < pow2(stage – 1); i++)
        Indx[pow2(stage – 1) + i] = Indx[i] + pow2(NumOfStage – stage);
 
  /* インデックスの並び替え */
  for (k = 0; k < N; k++)
    if (Indx[k] > k)
    {
      real = x_real[Indx[k]];
      imag = x_imag[Indx[k]];
      x_real[Indx[k]] = x_real[k];
      x_imag[Indx[k]] = x_imag[k];
      x_real[k] = real;
      x_imag[k] = imag;
    }
 
  free(Indx);
}

static void IFFT(double x_real[], double x_imag[], int N)
{
  int i;
  int j;
  int k;
  int n;
  int m;
  int r;
  int stage;
  int NumOfStage;
  int *Indx;
  double a_real;
  double a_imag;
  double b_real;
  double b_imag;
  double c_real;
  double c_imag;
  double real;
  double imag;
 
  /* IFFTの段数 */
  NumOfStage = log2(N);
 
  /* バタフライ計算 */
  for (stage = 1; stage <= NumOfStage; stage++)
    for (i = 0; i < pow2(stage – 1); i++)
      for (j = 0; j < pow2(NumOfStage – stage); j++)
      {
        n = pow2(NumOfStage – stage + 1) * i + j;
        m = pow2(NumOfStage – stage) + n;
        r = pow2(stage – 1) * j;
        a_real = x_real[n];
        a_imag = x_imag[n];
        b_real = x_real[m];
        b_imag = x_imag[m];
        c_real = cos((2.0 * M_PI * r) / N);
        c_imag = sin((2.0 * M_PI * r) / N);
        if (stage < NumOfStage)
        {
          x_real[n] = a_real + b_real;
          x_imag[n] = a_imag + b_imag;
          x_real[m] = (a_real – b_real) * c_real – (a_imag – b_imag) * c_imag;
          x_imag[m] = (a_imag – b_imag) * c_real + (a_real – b_real) * c_imag;
        }
        else
        {
          x_real[n] = a_real + b_real;
          x_imag[n] = a_imag + b_imag;
          x_real[m] = a_real – b_real;
          x_imag[m] = a_imag – b_imag;
        }
      }
 
  /* インデックスの並び替えのためのテーブルの作成 */
  Indx = calloc(N, sizeof(int));
  for (stage = 1; stage <= NumOfStage; stage++)
      for (i = 0; i < pow2(stage – 1); i++)
          Indx[pow2(stage – 1) + i] = Indx[i] + pow2(NumOfStage – stage);
 
  /* インデックスの並び替え */
  for (k = 0; k < N; k++)
    if (Indx[k] > k)
    {
      real = x_real[Indx[k]];
      imag = x_imag[Indx[k]];
      x_real[Indx[k]] = x_real[k];
      x_imag[Indx[k]] = x_imag[k];
      x_real[k] = real;
      x_imag[k] = imag;
    }
 
  /* 計算結果をNで割る */
  for (k = 0; k < N; k++)
  {
    x_real[k] /= N;
    x_imag[k] /= N;
  }
 
  free(Indx);
}

static int log2(int x)
{
  int y;

  for(y=0 ; x>1 ; y++)    x >>= 1;
  return y;
}

static int pow2(int x)
{
  int y;

  if (x == 0)   y = 1;
  else            y = 2 << (x – 1);
  return y;
}

static double max(const double *TMData, int N)
{
    int i;
    double MaxVl=0;

    for(i=0 ; i<N ; i++)
        if(MaxVl<TMData[i])
            MaxVl=TMData[i];

    return MaxVl;
}

————————————————————–

カテゴリー: 信号処理 パーマリンク