上课和平时实验是用的FFT都是matlab自带的,最近需要在一个嵌入式处理器中对采集的信号进行FFT实时处理,所以想先在matlab中自己写一个fft,然后在底层去实现fft算法。
下面是matlab的代码(使用函数进行了封装):
function [ ret_val ] = FFT1024( vector )
%UNTITLED8 此处显示有关此函数的摘要
% 此处显示详细说明
%UNTITLED7 此处显示有关此函数的摘要
% 此处显示详细说明
%======================================
%ret_val 为fft变换后返回的频域序列
%N 为点数
%vector 为变换前的序列
%======================================
vector_size = size(vector);
N = vector_size(2);
c = zeros(1,N);
%
%变址运算
%
j1 = 0;
for i = 1 : N
if i < j1 + 1
tmp = vector(j1 + 1);
vector(j1 + 1) = vector(i);
vector(i) =tmp;
end
k = N / 2;
while k <= j1
j1 = j1 - k;
k = k / 2;
end
j1 = j1 + k;
end
%
%蝶形运算
%
%%%%%%%计算 N 的
dig = 0;
k = N;
while k > 1
dig = dig + 1;
k = k / 2;
end
%%%%%%
% m 为级; dist 为蝶形运两点的距离; n 为蝶形运算组数
%
n = N / 2;
for m = 1 : dig
dist = 2 ^ (m - 1);
idx = 1;
for i = 1 : n
idx1 = idx;
for j1 = 1 : N / (2 * n)
r = (idx - 1) * 2 ^ (dig - m);
coef = exp(j * (-2 * pi * r / N));
tmp = vector(idx);
vector(idx) = tmp + vector(idx + dist) * coef;
vector(idx + dist) = tmp - vector(idx + dist) * coef;
idx = idx + 1;
end
idx = idx1 + 2 * dist;
end
n = n / 2;
end
ret_val = vector;
end
测试代码:
clc;
clear;
AD_fre = 1023 %修改采样频率和采样点数,采样时间没有修改,通过修改改值可以观察周期截取对信号频谱的影响
Adc=2; %直流分量幅度
A1=3; %频率F1信号的幅度
A2=1.5; %频率F2信号的幅度
F1=50; %信号1频率(Hz)
F2=75; %信号2频率(Hz)
Fs=AD_fre; %采样频率(Hz)
P1=0; %信号1相位(度)
P2=0; %信号相位(度)
N=AD_fre; %采样点数
t=[0:1/Fs:N/Fs]; %采样时刻
%信号
vector=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
%显示原始信号
figure(1)
plot(vector);
title('原始信号');
Y = FFT1024(vector); %做FFT变换
Ayy = (abs(Y)); %取模
Ayy=Ayy/(N/2); %换算成实际的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)*Fs/N; %换算成实际的频率值,Fn=(n-1)*Fs/N
figure(2)
stem(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果
title('幅度-频率曲线图');
%
% Pyy=[1:N/2];
% for i=1:N/2
% Pyy(i)=angle(Y(i)); %计算相位
% Pyy(i)=Pyy(i)*180/pi; %换算为角度
% end;
% subplot(414);stem(F(1:N/2),Pyy(1:N/2)); %显示相位图
% title('相位-频率曲线图');