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

オススメの参考書

C言語学習中だけど分からないことが多くて挫折しそう...という方には、下記の「スッキリわかるC言語入門」がオススメです!

まず学習を進める上で、参考書は2冊持っておくことをオススメします。この理由は下記の2つです。

  • 参考書によって、解説の仕方は異なる
  • 読み手によって、理解しやすい解説の仕方は異なる

ある人の説明聞いても理解できなかったけど、他の人からちょっと違った観点での説明を聞いて「あー、そういうことね!」って簡単に理解できた経験をお持ちの方も多いのではないでしょうか?

それと同じで、1冊の参考書を読んで理解できない事も、他の参考書とは異なる内容の解説を読むことで理解できる可能性があります。

なので、参考書は2冊持っておいた方が学習時に挫折しにくいというのが私の考えです。

特に上記の「スッキリわかるC言語入門」は、他の参考書とは違った切り口での解説が豊富で、他の参考書で理解できなかった内容に対して違った観点での解説を読むことができ、オススメです。題名の通り「なぜそうなるのか?」がスッキリ理解できるような解説内容にもなっており、C言語入門書としてもかなり分かりやすい参考書だと思います。

もちろんネット等でも色んな観点からの解説を読むことが出来ますので、分からない点は別の人・別の参考書の解説を読んで解決していきましょう!もちろん私のサイトも参考にしていただけると嬉しいです!

入門用のオススメ参考書は下記ページでも紹介していますので、こちらも是非参考にしていただければと思います。

https://daeudaeu.com/c_reference_book/

同じカテゴリのページ一覧を表示