matlab如何用fft

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)|')

运行结果。

% 基于Matlab的时间抽取基2FFT算法

function y=myditfft(x)

%本程序对输纤缓拆森入序列实现DIT-FFT基2算法，点数取大于等于长度的2的幂次

%------------------------------------

%Leo's fft program(改编网上的一个程序)

%------------------------------------

m=log2(2^nextpow2(length(x))) %求的x长度对应的2的最低毁御模幂次m

N=2^m

if length(x)<N

x=[x,zeros(1,N-length(x))] %若长度不是2的幂，补0到2的整数幂

end

x

%--------------------------------------------------------------------------

%对输入序列进行倒序

%如果输入序列的自然顺序号I用二进制数（例如n2n1n0）表示

%则其倒位序J对应的二进制数就是（n0n1n2），这样，在原来自然顺序时应该放x(I)的

%单元，现在倒位序后应放x(J)。

%--------------------------------------------------------------------------

%以下程序相当于以下程序:

%nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1 %求1：2^m数列的倒序

%y=x(nxd)%将倒序排列作为初始值

%--------------------------------------------------------------------------

NV2=N/2

NM1=N-1

I=0

J=0

while I<NM1

if I<J

T=x(J+1)

x(J+1)=x(I+1)

x(I+1)=T

end

K=NV2

while K<=J

J=J-K

K=K/2

end

J=J+K

I=I+1

end

x

%--------------------------------------------------------------------------

%以下程序解释:

%第一级从x(0)开始,跨接一阶蝶形,再取每条对称

%第二级从x(0)开始,跨接两阶蝶形,再取每条对称

%第m级从x(0)开始,跨接2^(m-1)阶蝶形,再取每条对称....

%--------------------------------------------------------------------------

for mm=1:m%将DFT做m次基2分解，从左到右，对每次分解作DFT运算

Nmr=2^mm

u=1 %旋转因子u初始化

WN=exp(-j*2*pi/Nmr) %本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)

for n=1:Nmr/2 %本次跨越间隔内的各次碟形运算

for k=n:Nmr:N %本次碟形运算的跨越间隔为Nmr=2^mm

kp=k+Nmr/2%确定碟形运算的对应单元下标(对称性)

t=x(kp)*u %碟形运算的乘积项

x(kp)=x(k)-t %碟形运算的加法项

x(k)=x(k)+t

end

u=u*WN%修改旋转因子，多乘一个基本DFT因子WN

end

end

y=x %输出

% 下面的程序里Pn 存的就是基波相位 如果求的是谐波相位，稍微修改即可

x = load('data.dat')%load 数据

fs=10000% 采样频率乎辩裤，自己根据实际情况设置

N=length(x)% x 是待分析的数据

n=1:N

%1-FFT

X=fft(x)% FFT

X=X(1:N/2)

Xabs=abs(X)

Xabs(1) = 0%直流分量岁简置0

[Amax,index]=max(Xabs)

if(Xabs(index-1) >Xabs(index+1))

a1 = Xabs(index-1) / Xabs(index)

r1 = 1/(1+a1)

k01 = index -1

else

a1 = Xabs(index) / Xabs(index+1)

r1 = 1/(1+a1)

k01 = index

end

Fn = (k01+r1-1)*fs/N%基灶配波频率

An = 2*pi*r1*Xabs(k01)/(N*sin(r1*pi))%基波幅值

Pn = phase(X(k01))-pi*r1%基波相位 单位弧度

Pn = mod(Pn(1),pi)

---