APP
【經(jīng)典算法實現(xiàn) 40】一維傅里葉變換 及 逆變換代碼實現(xiàn)
本文鏈接：《??【經(jīng)典算法實現(xiàn) 40】一維傅里葉變換 及 逆變換代碼實現(xiàn)??》
一、一維傅里葉變換 及 其逆變換
一維傅里葉變換：
F ( u ) = ∑ M = 0 M ? 1 f ( x ) e ? j u x 2 π M ， u = 0 , 1 , 2 , ? ? , M ? 1 F(u) = sum_{M=0}^{M-1} f(x) e^{-j u x frac{2pi} M}， u=0,1,2,cdots,M-1 F(u)=M=0∑M?1f(x)e?juxM2π，u=0,1,2,?,M?1
一維傅里葉逆變換：
f ( x ) = 1 M ∑ u = 0 M ? 1 F ( u ) e j u x 2 π M ， x = 0 , 1 , 2 , ? ? , M ? 1 f(x) = {frac 1 M} sum_{u=0}^{M-1} F(u) e^{j u x frac{2pi} M}， x=0,1,2,cdots,M-1 f(x)=M1?u=0∑M?1?F(u)ejuxM2π?，x=0,1,2,?,M?1
代入歐拉公式:
e j x = cos ? ( x ) + j sin ? ( x ) e^{jx} = cos(x) + j sin(x) ejx=cos(x)+jsin(x)
e ? j x = cos ? ( x ) ? j sin ? ( x ) e^{-jx} = cos(x) - j sin(x) e?jx=cos(x)?jsin(x)
e j π + 1 = 0 e^{jpi} + 1 = 0 ejπ+1=0
sin ? ( x ) = e i x ? e ? j x 2 i sin(x) = frac {e^{ix} - e^{-jx}} {2i} sin(x)=2ieix?e?jx?
cos ? ( x ) = e i x + e ? j x 2 i cos(x) = frac {e^{ix} + e^{-jx}} {2i} cos(x)=2ieix+e?jx?
我們可以得到如下：
一維傅里葉變換：
F ( u ) = ∑ M = 0 M ? 1 f ( x ) ( e ? j u x 2 π M ) ， u = 0 , 1 , 2 , ? ? , M ? 1 ? F ( u ) = ∑ M = 0 M ? 1 f ( x ) ( cos ? ( u x 2 π M ) ? j sin ? ( u x 2 π M ) ) ， u = 0 , 1 , 2 , ? ? , M ? 1 F(u) = sum_{M=0}^{M-1} f(x) left( e^{-j u x frac{2pi} M} ight)， u=0,1,2,cdots,M-1 \ Rightarrow F(u) = sum_{M=0}^{M-1} f(x) left( cos(u x frac{2pi} M) - j sin(u x frac{2pi} M) ight)， u=0,1,2,cdots,M-1 F(u)=M=0∑M?1?f(x)(e?juxM2π?)，u=0,1,2,?,M?1?F(u)=M=0∑M?1?f(x)(cos(uxM2π?)?jsin(uxM2π?))，u=0,1,2,?,M?1
一維傅里葉逆變換：
f ( x ) = 1 M ∑ u = 0 M ? 1 F ( u ) e j u x 2 π M ， x = 0 , 1 , 2 , ? ? , M ? 1 ? f ( x ) = 1 M ∑ u = 0 M ? 1 F ( u ) ( cos ? ( u x 2 π M ) + j sin ? ( u x 2 π M ) ) ， x = 0 , 1 , 2 , ? ? , M ? 1 f(x) = {frac 1 M} sum_{u=0}^{M-1} F(u) e^{j u x frac{2pi} M}， x=0,1,2,cdots,M-1 \ Rightarrow f(x) = {frac 1 M} sum_{u=0}^{M-1} F(u) left( cos(u x frac{2pi} M) + j sin(u x frac{2pi} M) ight) ， x=0,1,2,cdots,M-1 f(x)=M1?u=0∑M?1?F(u)ejuxM2π?，x=0,1,2,?,M?1?f(x)=M1?u=0∑M?1?F(u)(cos(uxM2π?)+jsin(uxM2π?))，x=0,1,2,?,M?1
在計算每個 F ( u ) F(u) F(u)時，
由于 一維傅里葉逆變換中的 u u u 是一個復數(shù)，
假設復數(shù) F ( u ) = a + b j F(u) = a + b j F(u)=a+bj，將 u u u 代入逆變換中，得到如下：
f ( x ) = 1 M ∑ u = 0 M ? 1 F ( u ) ( cos ? ( u x 2 π M ) + j sin ? ( u x 2 π M ) ) ， x = 0 , 1 , 2 , ? ? , M ? 1 f(x) = {frac 1 M} sum_{u=0}^{M-1} F(u) left( cos(u x frac{2pi} M) + j sin(u x frac{2pi} M) ight) ， x=0,1,2,cdots,M-1 f(x)=M1?u=0∑M?1?F(u)(cos(uxM2π?)+jsin(uxM2π?))，x=0,1,2,?,M?1
對于每個 F ( u ) F(u) F(u)，實部虛部計算結果如下：
f ( x ) = ( a + b j ) ( cos ? ( u x 2 π M ) + j sin ? ( u x 2 π M ) ) f(x) = ( a + b j) left( cos(u x frac{2pi} M) + j sin(u x frac{2pi} M) ight) f(x)=(a+bj)(cos(uxM2π?)+jsin(uxM2π?))
f ( x ) = ( a cos ? ( u x 2 π M ) + b j cos ? ( u x 2 π M ) + a j sin ? ( u x 2 π M ) + b j ? j sin ? ( u x 2 π M ) ) f(x) =left( acos(u x frac{2pi} M) + bjcos(u x frac{2pi} M) + aj sin(u x frac{2pi} M) + bjcdot j sin(u x frac{2pi} M) ight) f(x)=(acos(uxM2π?)+bjcos(uxM2π?)+ajsin(uxM2π?)+bj?jsin(uxM2π?))
分離虛部，實部結果如下：
f ( x ) . r e a l = ( a cos ? ( u x 2 π M ) ? b sin ? ( u x 2 π M ) ) f(x).real =left( acos(u x frac{2pi} M) - b sin(u x frac{2pi} M) ight) f(x).real=(acos(uxM2π?)?bsin(uxM2π?))
f ( x ) . i m a g = ( a sin ? ( u x 2 π M ) + b cos ? ( u x 2 π M ) ) ? j f(x).imag =left( a sin(u x frac{2pi} M) + bcos(u x frac{2pi} M) ight) cdot j f(x).imag=(asin(uxM2π?)+bcos(uxM2π?))?j
二、C代碼實現(xiàn)
// 一維傅 里葉變換代碼實現(xiàn)
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define NUM 15
typedef struct Complex_{
double x;
double y;
}Complex_;
float f[NUM];
Complex_ F[NUM];
Complex_ IF[NUM];
// 一維傅里葉變換
void DFT(void)
{
int i=0, j=0;
double tmp;
clock_t start, end;
start = clock();
srand(time(NULL));
// 初始化 F[NUM]
for(i = 0; i<NUM; i++){
//f[i] = sin(M_PI * i/(NUM*2));
f[i] = rand()%99;
F[i].x = 0; // 實數(shù)
F[i].y = 0; // 虛數(shù)
}
// M 求和
for(j = 0; j<NUM; j++)
{
// 計算每個數(shù)的 傅里葉變換結果
for(i = 0; i<NUM; i++)
{
tmp = -2 * M_PI * j * i / NUM;
F[j].x += f[i] * cos(tmp); // 實數(shù)
F[j].y += f[i] * sin(tmp); // 虛數(shù)
//printf("f[%d]=%f a[%d]=%9f b[%d]=%9f ",i, f[i], j, F[j].x, j, F[j].y);
}
}
end = clock();
printf(" 一維傅里葉變換 完畢, 使用時間:%lf NUM:%d ", (double)(end-start)/CLOCKS_PER_SEC, NUM);
for(i = 0; i<NUM; i++){
printf("f[%d]=%9f ---> F[%d]=%12f + %12fj ", i, f[i], i, F[i].x, F[i].y);
}
}
// 一維傅里葉逆變換
void IDFT(void)
{
int i=0, j=0;
double tmp;
clock_t start, end;
start = clock();
// 初始化 IF[NUM]
for(i = 0; i<NUM; i++){
IF[i].x = 0; // 實數(shù)
IF[i].y = 0; // 虛數(shù)
}
// M 求和
for(j = 0; j<NUM; j++)
{
// 計算每個數(shù)的 傅里葉變換結果
for(i = 0; i<NUM; i++)
{
tmp = 2 * M_PI * j * i / NUM;
IF[j].x += F[i].x * cos(tmp) - F[i].y * sin(tmp); // 實數(shù)
IF[j].y += F[i].x * sin(tmp) + F[i].y * cos(tmp); // 虛數(shù)
//printf("F[%d].x=%9f F[%d].y=%9f ---> IF[%d].x=%9f IF[%d].y=%9f ", i, F[i].x, i, F[i].y, j, IF[j].x, j, IF[j].y );
}
IF[j].x /= NUM;
IF[j].y /= NUM;
//printf("F[%d].x=%9f F[%d].y=%9f ---> IF[%d].x=%9f IF[%d].y=%9f ", j, F[j].x, j, F[j].y, j, IF[j].x, j, IF[j].y );
}
end = clock();
printf(" 一維傅里葉逆變換 完畢, 使用時間:%lf NUM:%d ", (double)(end-start)/CLOCKS_PER_SEC, NUM);
for(i = 0; i<NUM; i++){
printf("IF[%d] = %9f + %9fj ", i, IF[i].x, IF[i].y);
}
}
int main(void)
{
DFT();
printf(" ");
IDFT();
printf(" ");
return 0;
}
三、奇數(shù)位的變換結果
一維傅里葉變換 完畢, 使用時間:0.000000 NUM:15
f[ 0]=61.000000 ---> F[ 0]= 879.000000 + 0.000000j
f[ 1]=82.000000 ---> F[ 1]= 90.893436 + 13.541617j
f[ 2]=65.000000 ---> F[ 2]= -79.685020 + 5.178500j
f[ 3]=48.000000 ---> F[ 3]= 105.981226 + -28.670453j
f[ 4]=49.000000 ---> F[ 4]= -112.202453 + -67.242233j
f[ 5]=96.000000 ---> F[ 5]= -6.000000 + 6.928203j
f[ 6]=48.000000 ---> F[ 6]= -40.481226 + -21.662298j
f[ 7]=26.000000 ---> F[ 7]= 59.494037 + -155.029158j
f[ 8]= 1.000000 ---> F[ 8]= 59.494037 + 155.029158j
f[ 9]=88.000000 ---> F[ 9]= -40.481226 + 21.662298j
f[10]=45.000000 ---> F[10]= -6.000000 + -6.928203j
f[11]=88.000000 ---> F[11]= -112.202453 + 67.242233j
f[12]=44.000000 ---> F[12]= 105.981226 + 28.670453j
f[13]=89.000000 ---> F[13]= -79.685020 + -5.178500j
f[14]=49.000000 ---> F[14]= 90.893436 + -13.541617j
一維傅里葉逆變換 完畢, 使用時間:0.000000 NUM:15
IF[ 0] = 61.000000 + 0.000000j
IF[ 1] = 82.000000 + -0.000000j
IF[ 2] = 65.000000 + -0.000000j
IF[ 3] = 48.000000 + -0.000000j
IF[ 4] = 49.000000 + -0.000000j
IF[ 5] = 96.000000 + 0.000000j
IF[ 6] = 48.000000 + -0.000000j
IF[ 7] = 26.000000 + 0.000000j
IF[ 8] = 1.000000 + 0.000000j
IF[ 9] = 88.000000 + -0.000000j
IF[10] = 45.000000 + -0.000000j
IF[11] = 88.000000 + 0.000000j
IF[12] = 44.000000 + 0.000000j
IF[13] = 89.000000 + -0.000000j
IF[14] = 49.000000 + -0.000000j
四、偶數(shù)位的變換結果
一維傅里葉變換 完畢, 使用時間:0.000000 NUM:16
f[ 0]=74.000000 ---> F[ 0]= 907.000000 + 0.000000j
f[ 1]=83.000000 ---> F[ 1]= 74.197457 + -62.069951j
f[ 2]=48.000000 ---> F[ 2]= 62.468037 + 55.887302j
f[ 3]=34.000000 ---> F[ 3]= 81.074579 + 0.703911j
f[ 4]=96.000000 ---> F[ 4]= 35.000000 + 64.000000j
f[ 5]=29.000000 ---> F[ 5]= -53.802501 + -55.656479j
f[ 6]=87.000000 ---> F[ 6]= -70.468037 + -118.112698j
f[ 7]=53.000000 ---> F[ 7]= -21.469535 + 121.569659j
f[ 8]=54.000000 ---> F[ 8]= 63.000000 + -0.000000j
f[ 9]=41.000000 ---> F[ 9]= -21.469535 + -121.569659j
f[10]=21.000000 ---> F[10]= -70.468037 + 118.112698j
f[11]=75.000000 ---> F[11]= -53.802501 + 55.656479j
f[12]=36.000000 ---> F[12]= 35.000000 + -64.000000j
f[13]=26.000000 ---> F[13]= 81.074579 + -0.703911j
f[14]=69.000000 ---> F[14]= 62.468037 + -55.887302j
f[15]=81.000000 ---> F[15]= 74.197457 + 62.069951j
一維傅里葉逆變換 完畢, 使用時間:0.000000 NUM:16
IF[ 0] = 74.000000 + 0.000000j
IF[ 1] = 83.000000 + -0.000000j
IF[ 2] = 48.000000 + -0.000000j
IF[ 3] = 34.000000 + -0.000000j
IF[ 4] = 96.000000 + -0.000000j
IF[ 5] = 29.000000 + -0.000000j
IF[ 6] = 87.000000 + 0.000000j
IF[ 7] = 53.000000 + -0.000000j
IF[ 8] = 54.000000 + -0.000000j
IF[ 9] = 41.000000 + -0.000000j
IF[10] = 21.000000 + -0.000000j
IF[11] = 75.000000 + -0.000000j
IF[12] = 36.000000 + 0.000000j
IF[13] = 26.000000 + 0.000000j
IF[14] = 69.000000 + 0.000000j
IF[15] = 81.000000 + -0.000000j
從上面的奇數(shù) 和 偶數(shù)位的傅里葉運算結果中，大家有沒有發(fā)現(xiàn)什么規(guī)律呢？
本文摘自 ：https://blog.51cto.com/c
