作者HuangJC (吹笛牧童)
看板Math
標題Re: [問題] 還有人記得傅立葉變換嗎？
時間Mon Oct 21 19:05:04 2019
直接寫程式來湊答案
列出主程式
#define M_PI 3.14159265358979323846
#define SAMPLE_RATE 4
#define FFT_N 4
complex<double> s[FFT_N] = {0};
void test()
{
int i;
int sampleRate = SAMPLE_RATE; //Hz
double deltaDeg = 2 * M_PI / sampleRate;
TRACE("input:\n");
for(i = 0; i < sampleRate; i++) {
double x = deltaDeg * i;
s[i] = 1 + sin(x);// + 2 * sin(2 * x) + 3 * sin(3 * x);
TRACE("s[%d]=%f\n", i, s[i]);
}
FFT(s, FFT_N);
TRACE("\nafter FFT\n");
for(i = 0; i < FFT_N; i++) {
double real = s[i].real();
double imag = s[i].imag();
TRACE("s[%d]=(%f,%f) = %f deg %f\n",
i, real, imag, sqrt(real * real + imag * imag),
atan(imag / real) / M_PI * 180 );
}
}
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 111.82.254.91 (臺灣)
※ 文章網址: https://webptt.com/m.aspx?n=bbs/Math/M.1571655907.A.7B7.html
以上程式 s[i] = 1 + sin(x);
因為用的是 Visual C++ 6.0 with MFC
所以 s[i] 這種寫法，是直接填入實部
我用的測試函數是含直流成份 1， 1HZ sin 波，相角 0 （超級單純）
在此我們以 4HZ 做取樣，填入 size 為 4 的陣列中
#define SAMPLE_RATE 4
#define FFT_N 4
來看輸出
input:
s[0]=1.000000
s[1]=2.000000
s[2]=1.000000
s[3]=0.000000
after FFT
s[0]=(4.000000,0.000000) = 4.000000 deg 0.000000
s[1]=(-0.000000,-2.000000) = 2.000000 deg 90.000000
s[2]=(0.000000,0.000000) = 0.000000 deg -1.#IND00
s[3]=(-0.000000,2.000000) = 2.000000 deg -90.000000
s[0] 是直流成份，只有實部，其值為 4
s[1] 是 1hz, 只有虛部，其值為 -2
這什麼意思？
https://www.itread01.com/content/1544320336.html
》在輸出訊號而言，相鄰兩個取樣點的間隔為頻率解析度fs/N，
》其中fs為取樣頻率，N為輸入訊號的取樣點數
所以頻率解析度 = 4 / 4 = 1
果然如此，是這樣解讀頻域的 XD
但有沒有文件能告訴我，為什麼 s[0] 要除以 4，才變成直流成份 1
而 s[1] 要除以 2，才變成振幅 1
而且為什麼它是表達成 (0, -2），這是 -90度喔！ 但我明明原函式是 0度
（程式運算 atan 時，在 90 附近會有誤差，差個負號，結果角度錯了）
※ 編輯: HuangJC (111.82.254.91 臺灣), 10/21/2019 19:28:48
我們改用 s[i] = 1 + sin(x + 30.0 / 180 * M_PI);
相角 30 度當成輸入波型去做好了
after FFT
s[0]=(4.000000,0.000000) = 4.000000 deg 0.000000
s[1]=(1.000000,-1.732051) = 2.000000 deg -60.000000
s[2]=(0.000000,0.000000) = 0.000000 deg -1.#IND00
s[3]=(1.000000,1.732051) = 2.000000 deg 60.000000
就可以看到 s[1] 輸出是 -60 度，也就是再加 90度就可以還原原波型
為什麼要這樣迂迴啊...
※ 編輯: HuangJC (111.82.254.91 臺灣), 10/21/2019 19:31:04
好，那我現在來驗算 取樣 500 點，因為必需補到 2^N
所以陣列大小是 512 個，可以修改我的程式變成這樣來測
#define SAMPLE_RATE 500
#define FFT_N 512
還是用單純的波型
s[i] = 1 + sin(x);
after FFT
s[0]=(500.000000,0.000000) = 500.000000 deg 0.000000
s[1]=(6.641283,-253.006292) = 253.093442 deg -88.496361
s[2]=(-13.022571,6.307806) = 14.469823 deg -25.844379
s[3]=(-12.571770,1.665885) = 12.681663 deg -7.548286
輸出結果，直流成份要 除以 500
那麼 s[1] 的值除以 250，大約是 1
不過我們計算到 253，角度也不再是 -90 度
換言之 500 補到 512 還是有影響的
這東西只是'可參考'，但補 0 不是不影響原波型的！
※ 編輯: HuangJC (111.82.254.91 臺灣), 10/21/2019 19:44:18
1F：→ HuangJC : 接下來更挑戰了 Orz，明天再找同學確認 10/22 01:31
2F：→ HuangJC : 這套網路上抓來的 FFT，基本上處理最多 512點 10/22 01:32
3F：→ HuangJC : 可是我收到的檔案快六萬點，這怎麼處理！ 10/22 01:32