FFTで周波数がうまく取り出せません

Question

下記のプログラムはsinwave.wav（2000Hz×1秒間のファイル）から音声データを配列に格納して
その配列をFFTにかけるというプログラムです。
2000Hzの場所でパワースペクトルが上がることを期待していたのですが
なぜか900Hzあたりでパワースペクトルが上がってしまいます。
サンプル数のあたりが間違っているんでしょうか？
お力添えお願いします。

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

typedef struct

{

    int fs; //サンプリング周波数

    int bits; //量子化bit数

    int L; //データ長

} WAV_PRM;

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

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

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

double* audio_read(WAV_PRM* prm, char* filename);

double* audiodatasize;

int main(void) {

    WAV_PRM prm_in;
    double* data;
    FILE* audio = fopen("sinwave.wav", "rb");
    FILE* txt = fopen("data.txt", "w");
    char fn[] = "sinwave.wav";

    int i, Fr = 4096;
    double dt = 0.00005;
    double sc;
    double* x, * y;
    x = (double*)malloc(sizeof(double) * N);
    y = (double*)malloc(sizeof(double) * N);


    double fq;
    int pw = 0;
    printf("STG 01\n");
    sc = 2.0 * PI * dt;
    for (i = 0; i < N; i++) {
        if (x != NULL && y != NULL) {
            y[i] = 0.0;
        }
    }
    x = audio_read(&prm_in, fn);    //配列xにsinwave.wavの波形データ（2000Hz）を代入
    printf("STG 02\n");
    //  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 (x != NULL && y != NULL) {
            pw = sqrt(x[i] * x[i] + y[i] * y[i]) * 10; // パワースペクトル
        }
        fq = i / (dt * Fr);   //  周波数
        //printf("%f = %d / (%f * %d)\n", fq, i, dt, N);
        fprintf(fp, "%fHz  %d\n", fq, pw);
    }
    printf("STG 03\n");
    fclose(fp);
    fclose(audio);
    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;
        }
    }
}



double* audio_read(WAV_PRM* prm, char* filename)

{

    //変数宣言

    FILE* fp;

    int n;

    double* data;

    char header_ID[4];

    long header_size;

    char header_type[4];

    char fmt_ID[4];

    long fmt_size;

    short fmt_format;

    short fmt_channel;

    long fmt_samples_per_sec;

    long fmt_bytes_per_sec;

    short fmt_block_size;

    short fmt_bits_per_sample;

    char data_ID[4];

    long data_size;

    short data_data;



    //wavファイルオープン

    fp = fopen(filename, "rb");



    //wavデータ読み込み

    fread(header_ID, 1, 4, fp);

    fread(&header_size, 4, 1, fp);

    fread(header_type, 1, 4, fp);

    fread(fmt_ID, 1, 4, fp);

    fread(&fmt_size, 4, 1, fp);

    fread(&fmt_format, 2, 1, fp);

    fread(&fmt_channel, 2, 1, fp);

    fread(&fmt_samples_per_sec, 4, 1, fp);

    fread(&fmt_bytes_per_sec, 4, 1, fp);

    fread(&fmt_block_size, 2, 1, fp);

    fread(&fmt_bits_per_sample, 2, 1, fp);

    fread(data_ID, 1, 4, fp);

    fread(&data_size, 4, 1, fp);



    //パラメータ代入

    prm->fs = fmt_samples_per_sec;

    prm->bits = fmt_bits_per_sample;

    prm->L = data_size / 2;



    //音声データ代入

    data = (double*)calloc(prm->L, sizeof(double));
    audiodatasize = (double*)malloc(sizeof(double) * prm->L);

    for (n = 0; n < prm->L; n++) {

        fread(&data_data, 2, 1, fp);

        if (data != NULL) {
            data[n] = (double)data_data / 32768.0;
        }

    }



    fclose(fp);

    return data;

}

出力されたデータがこちらです。

～～～～～～
.
.
.
888.671875Hz  0
893.554688Hz  0
898.437500Hz  0
903.320313Hz  1
908.203125Hz  4
913.085938Hz  0
917.968750Hz  0
922.851563Hz  0
.
.
.
～～～～～～

##開発環境
Visual Studio 2019 C++

Comments

[退会済みユーザー]

読んでなさそうなので、もう一度貼っておきます。
[質問するときのヒント] https://teratail.com/help/question-tips
[推奨していない質問] https://teratail.com/help/avoid-asking
一通り読んでみてから、自分の質問がおかしくないかよく考えてみてください。

[jbpb0]

> 2000Hz×1秒間のファイル

> 2000Hzの場所でパワースペクトルが上がることを期待

2000Hzでサンプリングしたデータからは、1000Hz未満のスペクトルしか分かりません
参考
http://kinokorori.hatenablog.com/entry/20161207/1481117395

だからと言って、それは900Hzが上がる理由にはならないので、他にも何か間違いがあるのだとは思いますけど

Answer 1

fqを2.2倍することで解決しました。

Comments

[thkana]

どういう信号を入力したのか、「事実」についての情報が全くないのでなんとも言えませんが、多分その対応は間違っているのではないかと思います。
なぜWAVファイルから取得したサンプリング周波数が計算に反映されていないのですか?

「周波数」の解釈も、どうでしょう...FFTの出力は、0～fs/2が0～N/2にマッピングされ、N/2～N-1まではその折返しになる、とかわかっているのでしょうか。