FFTで20000Hzまで計算できるようにしたい

Question

FFTで20000Hzまで計算できるようにしたい

C++でFFTのプログラムを下記のサイトを参考にして作っています。
http://www-in.aut.ac.jp/~minemura/pub/Csimu/C/FFT.html
周波数は25000Hzくらいまで表示できるのですが
パワースペクトルが周波数と対応しません。
どなたかご教示お願い致します。

該当のソースコード

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#pragma warning(disable:4996)

FILE* fp;              // ファイル入出力

#define   PI  4.0*atan(1.0)  // 円周率

void S_fft(double* x, double* y, int, int);
#define N 32768  // N=データ総数(2,4,8,16,32,64,128,256,512,1024 のどれか)

int main(void) {
    int i, count, Fr = 128;
    double dt = 0.01;
    double sc;
    double* x, * y;
    x = (double*)malloc(sizeof(double) * N);
    y = (double*)malloc(sizeof(double) * N);


    double fq, pw;

    sc = 2.0 * PI * dt;
    for (i = 0; i < N; i++) {     // 1Hzに60Hzのノイズが重なった波形の例
        x[i] = sin(sc * i);
        y[i] = 0.0;
    }
    //  Fourier変換    
    S_fft(x, y, N, -1);

    //for (i = 0; i < Fr; i++) printf(" i=%d %f %f \n", i, x[i], y[i]);

    //    計算結果をファイルに格納
    fp = fopen("fft.txt", "w");     // ファイル名は fft.dat
    for (i = 0; i < N; i++) {
        if (i % 256 == 0) {
            pw = sqrt(x[i] * x[i] + y[i] * y[i]) * 100; // パワースペクトル
        }
        fq = i / (dt * Fr);   //  周波数
        //printf("%f = %d / (%f * %d)\n", fq, i, dt, N);
        fprintf(fp, "%fHz  %f \n", fq, pw);
    }
    fclose(fp);
    return 0;
}

/*********  Fourier変換および逆Fourier変換を行う副関数  ***************/
/*   入力:  x[n]: 被解析波形（振幅）  n; 離散データの数               */
/*          ff:   Fourier変換の場合は -1, 逆Fourier変換の場合は 1     */
/*   出力： ak[n]:cos関数の係数、 bk[n]:sin関数の係数                 */
/**********************************************************************/
void S_fft(double ak[], double bk[], int n, int ff) {
    int i, j, k, k1, num, nhalf, phi, phi0;
    static int rot[N];
    double s, sc, c, a0, b0, tmp;
    for (i = 0; i < n; i++) rot[i] = 0;

    nhalf = n / 2; num = n / 2; sc = 2.0 * PI / n;
    while (num >= 1) {
        for (j = 0; j < n; j += 2 * num) {
            phi = rot[j] / 2;    phi0 = phi + nhalf;
            c = cos(sc * phi);   s = sin(sc * phi * ff);
            for (k = j; k < j + num; k++) {
                k1 = k + num;
                a0 = ak[k1] * c - bk[k1] * s; b0 = ak[k1] * s + bk[k1] * c;
                ak[k1] = ak[k] - a0;      bk[k1] = bk[k] - b0;
                ak[k] = ak[k] + a0;       bk[k] = bk[k] + b0;
                rot[k] = phi;           rot[k1] = phi0;
            }
        }
        num = num / 2;
    }
    if (ff < 0) {
        for (i = 0; i < n; i++) {
            ak[i] /= n;  bk[i] /= n;
        }
    }

    for (i = 0; i < n - 1; i++) {
        if ((j = rot[i]) > i) {
            tmp = ak[i]; ak[i] = ak[j]; ak[j] = tmp;
            tmp = bk[i]; bk[i] = bk[j]; bk[j] = tmp;
        }
    }
}

試したこと

mallocを使ってみましたが、下記の警告が出て異常終了してしまいます。

警告

NULLポインター'x'を逆参照しています。
NULLポインター'y'を逆参照しています。

補足情報（FW/ツールのバージョンなど）

Visual Studio 2019

Comments

[1T2R3M4]

x[1024], y[1024]
rot[1024]
これらも変更してますか。

[Anfaenger]

変更してみました！

[thkana]

(質問が編集されましたが、以前のエラーは出なくなったのですか? ソースコードはなにも変更していないようですが、なにが変わったのですか?-最新の更新しか見ていなかった勘違いなのでここまで削除)

一致っていうのは、二つ以上のなにかとなにかが同じであることをいいます。
一つはこのプログラムの結果として、それとなにが一致しないと言っているのですか?

直流、正弦波、矩形波、鋸歯状波など、既知の波形を入力したらその結果は妥当なのですか?

[cateye]

>#define N 32768
・・・・
> double sc, x[N], y[N]; ←これ、スタック壊してませんか？

[Anfaenger]

>>thkanaさん
言葉足らずですみません。
周波数とそれに対応するパワースペクトルが一致しないという意味です。
Nをいじってサンプリング周波数？は変えられたと思うのですが、パワースペクトルの範囲が0〜99.9Hzになっているようで……

[Anfaenger]

>>cateyeさん
壊れてはいないと思いますが……
なぜですか？

[thkana]

私が、あなたがやったのと同じ比較をして「一致していない」ことを確認しようとするなら、
どうすれば比較する二つのデータを入手できますか?

[thkana]

もう一つ、比較対照しているデータが間違っていない、ということは確かなのですか?

[cateye]

>なぜですか？
遅くなりましたが・・・以下を参考に
https://teratail.com/questions/255528

[Anfaenger]

>>thkanaさん
データは以下のようになります。

0.000000Hz 0.001109
0.781250Hz 0.001109
1.562500Hz 0.001109
2.343750Hz 0.001109
3.125000Hz 0.001109
.
.
.

0.781250Hzのところで高いパワースペクトルが取り出せることを期待していたのですが
0になってしまいます。多分データがずれてるんだと思います。

追記：もしかするとcateyeさんの仰るようにスタックオーバーフロー？を起こしてるのかもしれません

[Anfaenger]

>>cateyeさん
mallocについても調べてみましたが初心者でよく分かりません……
お力添えをお願いします

Answer 1

FFTの周波数ってサンプリング周波数で正規化されてますから、単純にHzで表される話じゃないです。Nを増やしても「周波数」は変わりません。分解能が変わるだけです。

で、ソース中のコメントにN=データ総数(2,4,8,16,32,64,128,256,512,1024 のどれか)と書いてあります。何も考えずにNを1025以上にするなら異常動作が起こっても当然、でしょう。

内容は検討していませんが、1024で取っている配列の要素数をN以上に増やせばなんとかなりそうな気はします。

---

うわぁ、編集していっぱい書いたのが消えた...

書き直す気力もないので、簡略版。

    double dt = 0.01;
//<略>
    sc = 2.0 * PI * dt;
    for (i = 0; i < N; i++) {     // 1Hzに60Hzのノイズが重なった波形の例 <- ウソのコメント
        x[i] = sin(sc * i);
        y[i] = 0.0; //本題に関係ない
    }

つまり、sin(2pi0.01*i)だから100サンプル周期の正弦波です。32768/100=327.68ですから、327とか328番目(最初はDC成分なのでさらに+1する)を見てみると、

328番目 19.771515
329番目 41.983550
330番目 10.170375

ちゃんと出ているようにおもいます。

なにか問題がありますか? ナントカHzってのは全く意味不明なので無視していますが。