作者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/cn.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