C语言编写一个一维傅里叶函数

#include<stdio.h>

#include <math.h>

class complex //定义一个类,实现复数的所有 *** 作

{

double Real,Image//实部与虚部

public:

complex(double r="0",double i="0"){Real=rImage=i}

double GetR(){return Real}//取出实部

double GetI(){return Image} //取出虚部

complex operator + (complex &) //复数加法

complex operator - (complex &) //复数减法

complex operator * (complex &) //复数乘法

void operator =(complex &) //复数 赋值

}

complex complex::operator + (complex &c) //复数加法

{

complex t

t.Real=Real+c.Real

t.Image=Image+c.Image

return t

}

complex complex::operator - (complex &c) //复数减法

{

complex t

t.Real=Real-c.Real

t.Image=Image-c.Image

return t

}

complex complex::operator * (complex &c) //复数乘法

{

complex t

t.Real=Real*c.Real-Image*c.Image

t.Image=Real*c.Image+Image*c.Real

return t

}

void complex::operator = (complex &c) //复数 赋值

{

Real=c.Real

Image=c.Image

}

void fft(complex a[],int length,int jishu) //实现fft的函数

{

const double PI="3".141592653589793

complex u,Wn,t

int i,j,k,m,kind,distance,other

double tmp

for(i=0i<lengthi++) //实现倒叙排列

{

k="i"

j=0

for(m=0m<jishum++)

{

j="j"*2+k%2

k/=2

}

if(i<j)

{

t="a"

a=a[j]

a[j]=t

}

}

for(m=1m<=jishum++) //第m级蝶形运算,总级数为jishu

{

kind = (int)pow(2,m-1) //第m级有2^(m-1)种蝶形运算

distance = 2*kind //同种蝶形结相邻距离为2^m

u=complex(1,0)//旋转因子初始值为 1

tmp=PI/kind

Wn=complex(cos(tmp),-sin(tmp))//旋转因子Wn

for(j=0j<kindj++) //每种蝶形运算的起始点为j,共有kind种

{

for(i=ji<lengthi+=distance) //同种蝶形运算

{

other=i+kind//蝶形运算的两个因子对应单元下标的距离为2^(m-1)

t=a[other]*u// 蝶形运算的乘积项

a[other]=a-t //蝶形运算

a=a+t //蝶形运算

}

u="u"*Wn//修改旋转因子,多乘一个基本DFT因子WN

}

}

}

void main(void)

{

double a,b

complex x[8]//此程序以8点序列测试

printf("8点序列:\n")

for(int i="0"i<8i++) //初始化并输出原始序列

{

x=complex(i,i+1)

printf("x(%d) = %lf + %lf i\n",i+1,x.GetR(),x.GetI())

}

fft(x,8,3) //调用fft函数

printf("fft变换的结果为:\n")

for(i=0i<8i++) //输出结果

printf("X(%d)= %lf + %lf i\n",i+1,x.GetR(),x.GetI())

}

matlab自带的fft函数是快速傅里叶变换函数。主要用于降噪处理，通过使用傅里叶变换求噪声中隐藏的信号的频率分量。

该函数使用方法：

方法一：

Y = fft(X) 用快速傅里叶变换 (FFT) 算法计算 X 的离散傅里叶变换 (DFT)。

如果 X 是向量，则 fft(X) 返回该向量的傅里叶变换。

如果 X 是矩阵，则 fft(X) 将 X 的各列视为向量，并返回每列的傅里叶变换。

如果 X 是一个多维数组，则 fft(X) 将沿大小不等于 1 的第一个数组维度的值视为向量，并返回每个向量的傅里叶变换。

方法二：

Y = fft(X,n) 返回 n 点 DFT。如果未指定任何值，则 Y 的大小与 X 相同。

如果 X 是向量且 X 的长度小于 n，则为 X 补上尾零以达到长度 n。

如果 X 是向量且 X 的长度大于 n，则对 X 进行截断以达到长度 n。

如果 X 是矩阵，则每列的处理与在向量情况下相同。

如果 X 为多维数组，则大小不等于 1 的第一个数组维度的处理与在向量情况下相同。

我们通过下例，来了解fft函数使用过程：

第一步、指定信号的参数，采样频率为 1 kHz，信号持续时间为 1.5 秒。

Fs=1000；%采样频率

T=1/Fs；%采样周期

L=1500；%信号长度

t=（0:L-1）*T；%时间向量

第二步、构造一个信号，其中包含幅值为 0.7 的 50 Hz 正弦量和幅值为 1 的 120 Hz 正弦量。

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t)

第三步、用均值为零、方差为 4 的白噪声扰乱该信号。

X = S + 2*randn(size(t))

第四步、在时域中绘制含噪信号。通过查看信号 X(t) 很难确定频率分量。

plot(1000*t(1:50),X(1:50))

title('Signal Corrupted with Zero-Mean Random Noise')

xlabel('t (milliseconds)')，ylabel('X(t)')

第五步、计算信号的傅里叶变换。

Y = fft(X)

第六步、计算双侧频谱 P2， 计算单侧频谱 P1。

P2 = abs(Y/L)

P1 = P2(1:L/2+1)

P1(2:end-1) = 2*P1(2:end-1)

第七步、定义频域 f 并绘制单侧幅值频谱 P1

f = Fs*(0:(L/2))/L

plot(f,P1)

title('Single-Sided Amplitude Spectrum of X(t)')

xlabel('f (Hz)')，ylabel('|P1(f)|')

运行结果。

---