2次元離散フーリエ変換(DFT)をプログラミングしてみたけど・・・重すぎる!

周波数解析関連のプログラミングやフィルタ処理なんかを紹介しようと思い、まず周波数変換するための2次元離散フーリエ変換を実装してみました。

ですが、プログラムが終了してくれない・・・。

どうも離散フーリエ変換の部分のループが多すぎるらしい。500×500の画像ですら処理終了してくれなかった・・・。

高速フーリエ変換のすごさを実感しました。高速フーリエ変換を早速実装しようと思ったのですが2次元のものを作ろうと思うとなかなか難しい・・。

いつか完成して記事にしたいと思います!!!

#include "myJpeg.h"
#include <math.h>

#define PI 3.14

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

  RAWDATA_t raw, out;
  int m, n, c;
  int u, v;
  double *rdata, *idata;
  double rsum, isum;
  unsigned int ulen, vlen;

  char outname[256];
  FILE *fo;

  if(argc != 2){
    printf("ファイル名を引数に指定してください\n");
    return -1;
  }

  if(jpegFileReadDecode(&raw, argv[1]) == -1){
    printf("jpegFileReadDecode error\n");
    return -1;
  }

  ulen = raw.width;
  vlen = raw.height;

  rdata = (double*)malloc(sizeof(double) * ulen * vlen * raw.ch);
  if(rdata == NULL){
    printf("malloc rdata error\n");
    freeRawData(&raw);
    return -1;
  }

  idata = (double*)malloc(sizeof(double) * ulen * vlen * raw.ch);
  if(idata == NULL){
    printf("malloc idata error\n");
    freeRawData(&raw);
    free(rdata);
    return -1;
  }


  for(v = 0; v < vlen; v++){
    for(u = 0; u < ulen; u++){
      for(c = 0; c < raw.ch; c++){
        rsum = 0;
        isum = 0;
        for(n = 0; n < raw.height; n++){
         for(m = 0; m < raw.width; m++){
            rsum += raw.data[raw.ch * (m + n * raw.width) + c] * cos(2 * PI * (u * m / (double)raw.width + v * n / (double)raw.height));
            isum -= raw.data[raw.ch * (m + n * raw.width) + c] * sin(2 * PI * (u * m / (double)raw.width + v * n / (double)raw.height));
          }
       }
        rdata[raw.ch * (u + v * ulen) + c] = rsum;
        idata[raw.ch * (u + v * ulen) + c] = isum;
      }
    }
  }
  /* ここまで画像処理 */

  out.data = (unsigned char*)malloc(sizeof(unsigned char) * ulen * vlen * raw.ch);
  if(out.data == NULL){
    printf("malloc data error\n");
    freeRawData(&raw);
    free(rdata);
    free(idata);
    return -1;
  }
  out.width = ulen;
  out.height = vlen;
  out.ch = raw.ch;

  /* パワースペクトルを計算しようとしているがここまで実行されない・・・ */
  for(v = 0; v < vlen; v++){
    for(u = 0; u < ulen; u++){
      for(c = 0; c < out.ch; c++){
        out.data[out.ch * (u + v * ulen) + c] =
          rdata[raw.ch * (u + v * ulen) + c] * rdata[raw.ch * (u + v * ulen) + c] +
          idata[raw.ch * (u + v * ulen) + c] * idata[raw.ch * (u + v * ulen) + c];
      }
    }
  }

  sprintf(outname, "%s", "dft.jpeg");

  if(jpegFileEncodeWrite(&out, outname) == -1){
    printf("jpegFileEncodeWrite error\n");
    freeRawData(&raw);
    free(rdata);
    free(idata);
    return -1;
  }

  freeRawData(&raw);
  free(rdata);
  free(idata);
  freeRawData(&out);

  return 0;
}

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です