離散フーリエ変換をしたいです．

前提・実現したいこと

離散フーリエ変換をしたいのですが，現在，分割数N1とサンプリング周波数によってできたx配列の値があっていないとうまくフーリエ逆変換ができません．
サンプリング周波数が100Hzのとき1[s]までしているので配列の個数は100個になります．一周期で最低4つ値をとれば三角波になると思うのですが，うまくいきません．どなたかお願いします．
また，矩形波も作り方がいまいちわからないので，もしよければご指摘のほどお願いします．

発生している問題・エラーメッセージ

分割数と元のデータ数が同じ出ないと逆フーリエ変換ができない．

エラーメッセージ

該当のソースコード

c++
1//ヘッダーファイルです
2#define _USE_MATH_DEFINES
3#include <iostream>
4#include <vector>
5#include <cmath>
6#include <complex>
7
8#define N 1
9#define N1 100
10
11class Wave {
12private:
13    /* data */
14    void Wave_Output(std::vector<double> dt);//csvファイル出力(波)
15    void Fourier_Out(std::vector<double> dt_re, std::vector<double> dt_im);
16    void Fourier_R_Out(std::vector<double>dt_re, std::vector<double> dt_im);
17    
18public:
19    std::vector<double> x; //波
20    std::vector<double> x_r;//フーリエ逆変換で生成した波
21    std::vector<double> dt_re;
22    std::vector<double> dt_im;
23    std::vector<double> f_dt;
24    void test();
25    /*!---------------------------------------------------------
26    @brief 正弦波信号x(t)を生成する関数
27
28    double *dt 型の配列をキーボードから入力する
29    @return dt[]に入力した値が返る．
30    @attention 特になし
31    @par 更新履歴
32    - 2018/04/22
33    - 基本的な機能の実装(by)
34
35    */
36    void Sin_Wave(double A, double B, double f);
37    /*!---------------------------------------------------------
38    @brief 合成波1x(t)を生成する関数
39
40    double *dt 型の配列をキーボードから入力する
41    @return dt[]に入力した値が返る．
42    @attention 特になし
43    @par 更新履歴
44    - 2018/04/22
45    - 基本的な機能の実装(by)
46
47    */
48    void Comp_Wave1(double A1, double B1, double A2, double B2, double f);
49    /*!---------------------------------------------------------
50    @brief 合成波2x(t)を生成する関数
51
52    double *dt 型の配列をキーボードから入力する
53    @return dt[]に入力した値が返る．
54    @attention 特になし
55    @par 更新履歴
56    - 2018/04/22
57    - 基本的な機能の実装(by )
58
59    */
60    void Comp_Wave2(double A1, double B1, double A2, double B2, double C, double f);
61    /*!---------------------------------------------------------
62    @brief 矩形波x(t)を生成する関数
63
64    double *dt 型の配列をキーボードから入力する
65    @return dt[]に入力した値が返る．
66    @attention 特になし
67    @par 更新履歴
68    - 2018/04/22
69    - 基本的な機能の実装(by A)
70
71    */
72    void Rect_Wave(double A, int inter,double f);
73    /*!---------------------------------------------------------
74    @brief 離散フーリエ変換をする関数
75
76    double *dt 型の配列をキーボードから入力する
77    @return dt[]に入力した値が返る．
78    @attention 特になし
79    @par 更新履歴
80    - 2018/04/22
81    - 基本的な機能の実装(b)
82
83    */
84    void Fourier(int N_s);
85    /*!---------------------------------------------------------
86    @brief 離散フーリエ逆変換をする関数
87
88    double *dt 型の配列をキーボードから入力する
89    @return dt[]に入力した値が返る．
90    @attention 特になし
91    @par 更新履歴
92    - 2018/04/22
93    - 基本的な機能の実装(by )
94
95    */
96    void Reverse_Fourier(int N_r);
97    /*!---------------------------------------------------------
98    @brief パワースペクトルを計算する関数
99
100    double *dt 型の配列をキーボードから入力する
101    @return dt[]に入力した値が返る．
102    @attention 特になし
103    @par 更新履歴
104    - 2018/04/22
105    - 基本的な機能の実装(byA)
106
107    */
108    void Power_Spectral(double f);
109};

c++
1//関数ファイルです．
2#include "****.h"
3
4void Wave::Wave_Output(std::vector<double> dt) {
5    FILE *fp;
6    if ((fopen_s(&fp, "wave.csv", "w")) != 0) {
7        fprintf(stderr, "ファイルのオープンに失敗");
8    }
9    else {
10        for (int i = 0; i < dt.size(); i++) {
11            fprintf(fp, "%f\n", dt[i]);
12        }
13        fclose(fp);
14    }
15}
16
17void Wave::Fourier_Out(std::vector<double> dt_re, std::vector<double> dt_im) {
18    FILE *fp;
19    if ((fopen_s(&fp, "Fourier.csv", "w")) != 0) {
20        fprintf(stderr, "ファイルのオープンに失敗");
21    }
22    else {
23        for (int i = 0; i < dt_re.size(); i++) {
24            fprintf(fp, "%d,%f,%f\n", i * N / N1, dt_re[i], dt_im[i]);
25        }
26        fclose(fp);
27    }
28}
29
30void Wave::Fourier_R_Out(std::vector<double>dt_re, std::vector<double> dt_im) {
31    FILE *fp;
32    if ((fopen_s(&fp, "Fourier_R.csv", "w")) != 0) {
33        fprintf(stderr, "ファイルのオープンに失敗");
34    }
35    else {
36        for (int i = 0; i < dt_re.size(); i++) {
37            fprintf(fp, "%d,%f,%f\n", i * N / N1, dt_re[i], dt_im[i]);
38        }
39        fclose(fp);
40    }
41}
42
43void Wave::Power_Spectral(double f) {
44    FILE *fp;
45    if ((fopen_s(&fp, "Power_Spe.csv", "w")) != 0) {
46        fprintf(stderr, "ファイルのオープンに失敗");
47    }
48    else {
49        for (int i = 0; i < dt_re.size(); i++) {
50            fprintf(fp, "%f,%f\n", f_dt[i], sqrt(dt_re[i] * dt_re[i] + dt_im[i] * dt_im[i]));
51        }
52        fclose(fp);
53    }
54}
55void Wave::test() {
56    std::cout << "test" << '\n';
57}
58
59//基本周波数は
60void Wave::Sin_Wave(double A, double B, double f) {
61    for (double t = 0; t < N; t += 1.0 / f) {
62        /* code */
63        x.push_back(A*sin(2 * M_PI * B * t));
64    }
65    Wave out;
66    out.Wave_Output(x);
67}
68
69void Wave::Comp_Wave1(double A1, double B1, double A2, double B2, double f) {
70    for (double t = 0; t < N; t += 1.0 / f) {
71        /* code */
72        x.push_back(A1*sin(B1*t * 2 * M_PI) + A2 * cos(B2*t * 2 * M_PI));
73    }
74    Wave out;
75    out.Wave_Output(x);
76}
77
78void Wave::Comp_Wave2(double A1, double B1, double A2, double B2, double C, double f) {
79    for (double t = 0; t < N; t += 1.0 / f) {
80        /* code */
81        x.push_back(C + A1 * sin(B1*t * 2 * M_PI) + A2 * cos(B2*t * 2 * M_PI));
82    }
83    Wave out;
84    out.Wave_Output(x);
85}
86
87void Wave::Rect_Wave(double A, int inter,double f) {
88    double temp = A;
89    for (double t = 0; t < 2*3; t+= 1/f) {
90        if (t == 3) {
91            temp = -1 * temp;
92        }
93        x.push_back(temp);
94    }
95    Wave out;
96    out.Wave_Output(x);
97}
98
99void Wave::Fourier(int N_s) {
100    //実部と虚部を分けてフーリエ変換
101    std::vector<double> re(N_s);
102    std::vector<double> im(N_s);
103    for (int n = 0; n < N_s; n++) {
104        re[n] = 0.0;
105        im[n] = 0.0;
106        for (int k = 0; k < N_s; k ++) {
107            re[n] += x[k] * cos(2 * M_PI*k*n / N_s) / N_s;
108            im[n] += -x[k] * sin(2 * M_PI*k*n / N_s) / N_s;
109        }
110        f_dt.push_back(2 * M_PI*n / N_s);
111    }
112    Wave out;
113    out.Fourier_Out(re,im);
114    dt_re = re;
115    dt_im = im;
116}
117
118void Wave::Reverse_Fourier(int N_r) {
119    std::vector<double> re(N_r);
120    std::vector<double> im(N_r);
121    int k;
122    for (int n = 0; n < N_r; n++) {
123        re[n] = 0.0;
124        im[n] = 0.0;
125        for (k = 0; k < N_r; k++) {
126            re[n] += (dt_re[k] * cos(2 * M_PI*k*n / N_r)) - (dt_im[k] * sin(2 * M_PI*k*n / N_r));
127            im[n] += (dt_re[k] * sin(2 * M_PI*k*n / N_r)) + (dt_im[k] * cos(2 * M_PI*k*n / N_r));
128        }
129    }
130    Wave out;
131    out.Fourier_R_Out(re, im);
132}

c++
1//main文です
2#include "****.h"
3#include <iostream>
4#include <vector>
5
6
7
8int main() {
9
10    typedef int size_t;
11    //std::vector<double> v = A_Wave(1,-1,100.0);
12    Wave W;
13    //W.Comp_Wave2(1.0, 2.0, 2.0, 10.0, 3.0,100);
14    W.Sin_Wave(1.0,2.0,100);
15    //W.Rect_Wave(1,M_PI,10);
16    W.Fourier(N1);
17    W.Reverse_Fourier(N1);
18    W.Power_Spectral(100);
19    system("pause");
20    return 0;
21}
22

試したこと

お願いします．

chankane

2018/04/26 10:09

すみません、お答えできないと思いますが教えてください。W.Reverse_Fourier(N1) においてN1を任意の数字に指定したいということですか？

muscle1100

2018/04/26 11:19

そうですね。元から離散的に取っているところから、さらにデータを取るっていう方法（そもそもこれが間違っているかもしれません）をとっているので！

chankane

2018/04/26 11:33

すみません、まだあまりよくわかっていません (´･ω･`)ｺﾞﾒﾝﾖ　＞元から離散的に取っているところから、さらにデータを取るっていう方法（そもそもこれが間違っているかもしれません）をとっているので！　この一文をもっと具体的に説明できませんか？

muscle1100

2018/04/26 11:43

離散フーリエ変換は離散的ではない式を一定ごとに区切って離散的にすると思うのですが、そもそもの値が1/fごとの値を取っているから離散的だなと思いまして（説明下手ですみません）

chankane

2018/04/26 12:02

もしかしたら私の解釈が間違っているかもしれませんが、ご指摘致します。質問者様は離散フーリエ変換を少し誤解していらっしゃるように思いました。離散フーリエ変換自体にはアナログ信号をデジタル化する機能はございません。（そもそもの音声データが既に離散的）。もしも質問者様が「そんなことしってるよヽ(ﾟДﾟ)ﾉ」とお思いになったのでしたら、大変失礼いたしました。

muscle1100

2018/04/26 12:10

大丈夫です！ご指摘ありがとうございます

行動規範の内容に同意します

回答1件

ベストアンサー

すみません。当方では質問者様のおかれている状況が完全に理解できていないので、参考になりそうなURLをのせるだけになります。
お役に立てずにすみません(´･ω･`)
Pythonで短時間フーリエ変換（STFT）と逆変換
[STFTとISTFTの秘密の関係
フーリエ変換](http://d.hatena.ne.jp/aqua-nora/20111109/1320808234)

投稿2018/04/26 13:02

chankane

総合スコア139