求用C语言实现FFT变换的程序（见下面）

1）结果讨论

一，如果对信号进行同样点数N的FFT变换，采样频率fs越高，则可以分析越高频的信号；与此同时，采样频率越低，对于低频信号的频谱分辨率则越好。

二，假设采样点不在正弦信号的波峰、波谷、以及0电压处，频谱则会产生泄露（leakage）。

三，对于同样的采样率fs，提高FFT的点数N，则可提高频谱的分辨率。

四，如果采样频率fs小于2倍信号频率2*fs（奈圭斯特定理），则频谱分析结果会出错。

五，对于（二）中泄露现象，可以通过在信号后面补零点解决。

2）程序及注解如下

%清除命令窗口及变量

clc

clear all

%输入f、N、T、是否补零（补几个零）

f=input('Input frequency of the signal: f\n')

N=input('Input number of pointsl: N\n')

T=input('Input sampling time: T\n')

flag=input('Add zero too sampling signal or not? yes=1 no=0\n')

if(flag)

ZeroNum=input('Input nmber of zeros\n')

else

ZeroNum=0

end

%生成信号，signal是原信号。signal为采样信号。

fs=1/T

t=0:0.00001:T*(N+ZeroNum-1)

signal=sin(2*pi*f*t)

t2=0:T:T*(N+ZeroNum-1)

signal2=sin(2*pi*f*t2)

if (flag)

signal2=[signal2 zeros(1, ZeroNum)]

end

%画出原信号及采样信号。

figure

subplot(2,1,1)

plot(t,signal)

xlabel('Time(s)')

ylabel('Amplitude(volt)')

title('Singnal')

hold on

subplot(2,1,1)

stem(t2,signal2,'r')

axis([0 T*(N+ZeroNum) -1 1])

%作FFT变换，计算其幅值，归一化处理，并画出频谱。

Y = fft(signal2,N)

Pyy = Y.* conj(Y)

Pyy=(Pyy/sum(Pyy))*2

f=0:fs/(N-1):fs/24

subplot(2,1,2)

bar(f,Pyy(1:N/2))

xlabel('Frequency(Hz)')

ylabel('Amplitude')

title('Frequency compnents of signal')

axis([0 fs/2 0 ceil(max(Pyy))])

grid on

祝你好运！

我可以帮助你，你先设置我最佳答案后，我百度Hii教你。

1.首先建立信号模型并采样

%谐波分析仿真数据生成 CreatData.m

clc

fs=3000f0=50

N=1024

n=1:N

t=(n-1)/fs

m=13

Am=[380.4 0 2.035 0 3.156 0 1.042 0 4.291 0 16.13 0 2.016 ]

PH=[0 0 60 0 135 0 157.5 0 0 0 60 0 0]

x=zeros(1,N)

for k=1 : m

x=x+Am(k)*cos(2*pi*f0*k*t+PH(k)/180*pi)

end

figure(1)

plot(x)

2.进行FFT谐波分析

%直接用FFT进行谐波计算 Harmonic_Analysis.m

N=length(x)

f=50

m=13

n=1:N

Fn=zeros(9,3)% 谐波频率

An=zeros(9,3)% 幅值

Pn=zeros(9,3)% 初相角

for i=1 : m

Fn(i,1) = iFn(i,2) = f * i

An(i,1) = iAn(i,2) = 0

Pn(i,1) = iPn(i,2) = 0

end

fsN=fs/N

f0=50

%1-FFT

X=fft(x)

X=X(1:N/2)

Xabs=abs(X)

for i= 1 : m

[Amax,index]=TriFind(Xabs,floor((i*f0-15)/fsN),ceil((i*f0+15)/fsN))

if(index==-1)

Fn(i,3) = 0

An(i,3) = 0

Pn(i,3) = 0

else

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(i,3) = (k01+r1-1)*fs/N

An(i,3) = 2*pi*r1*Xabs(k01)/(N*sin(r1*pi))

Pn(i,3) = phase(X(k01))-pi*r1

Pn(i,3) = Pn(i,3)*180/pi

Pn(i,3) = mod(Pn(i,3),180)

end

end

Fn

An

Pn

figure(2)

stem(Xabs(1:N/2)*2/N)

dlmwrite('An.txt',An,'delimiter', '\t','precision', '%.10f','newline', 'pc')

dlmwrite('Fn.txt',Fn,'delimiter', '\t','precision', '%.10f','newline', 'pc')

dlmwrite('Pn.txt',Pn,'delimiter', '\t','precision', '%.10f','newline', 'pc')

3.完全手打的，先运行CreatData.m，再运行Harmonic_Analysis.m，应该可行，有问题再补充

4.少了个函数文件

新建个.m,命名为 TriFind，将以下内容粘贴

function [ Amax,index ] = TriFind(A,first,last)

%TRIFIND find x[k] if x[k-1]<x[k] and x[k]>x[k+1]

% return x[k] and the k

Amax=0

index=-1

for k=first+1:last-1

if A(k-1)<A(k) &&A(k)>A(k+1)

Amax=A(k)

index=k

end

end