
n年前做的了基2的
#include "stdafxh"
#include "complexh"
void main()
{
const double pi = 314159265;
int i; //循环变量
//////////////////////////////////////////////////////////////////
//信号生成:正弦信号,固定频率。生成1024点
const N = 64; //信号点数
const n =N/2; //半点数,用于循环控制
float S[N]; //信号数组
float temp;
temp = (float) (8 pi) / N; //出现两个周期
for(i=0; i<N; i++) S[i] = (float) sin(temp i);
//////////////////////////////////////////////////////////////////
//旋转因子数组生成
complex W[n];
temp = (float) (2 pi / N);
for(i=0; i<n; i++) {
W[i]SetReal((float)cos(temp i));
W[i]SetImag((float)-sin(temp i));
// cout<<W[i]<<endl;
}
///////////////////////////////////////////////////////////////////
//变换:基2-FFT
//先做蝶形运算,再对信号排序
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
//读入排序的信号
complex x;
x = new complex[Nsizeof(complex)];//计算用数组
///////////////////////////////////////////////////////////////////
//读入信号
for(i=0; i<N; i++) x[i]SetReal(S[i]);
///////////////////////////////////////////////////////////////////
//蝶形运算
int j, k; //循环变量
complex t;
for(i=n, k=1; i>=1; i>>=1, k<<=1) {
k %= n; //旋转因子调整
for(j=0; j<N;) {
t = x[j+i]; //蝶形运算
x[j+i] = x[j] - t;
x[j] = x[j] + t;
x[j+i] = x[j+i] W[(j%i)k];
j++;
if((j%i) == 0) j += i; //交叉
}
}
/////////////////////////////////////////////////////////////////////
//生成(以8点为例)形如0,4,2,6,1,5,3,7的排序数组
int y;
y = new int[Nsizeof(int)]; //调整计算用数组的排序用数组
for(i=0; i<N; i++) y[i] = 0;
for(i=1, k=n; i<=n; i<<=1, k>>=1)
for(j=0; j<2i; j++) y[i+j] = y[j] + k;
//////////////////////////////////////////////////////////////////////
//排序
complex X;
X = new complex[Nsizeof(complex)];
for(i=0; i<N; i++) X[y[i]] = x[i];
delete [] x; //销毁运算数组
delete [] y; //销毁排序数组
////////////////////////////////////////////////////////////////////////
//归一化
float a[N];
for(i=0; i<N; i++) a[i] = X[i]mod();
delete [] X;
temp = a[0];
for(i=1; i<N; i++) if(a[i]>temp) temp = a[i];
for(i=0; i<N; i++) a[i] =(float) a[i] / temp;
for(i=0; i<N; i+=4) cout<<a[i]<<" "<<a[i+1]<<" "<<a[i+2]<<" "<<a[i+3]<<endl;
}
在图象处理的广泛应用领域中,傅立叶变换起着非常重要的作用,具体表现在包括图象分析、图象增强及图象压缩等方面。
fftshift的作用正是让正半轴部分和负半轴部分的图像分别关于各自的中心对称。因为直接用fft得出的数据与频率不是对应的,fftshift可以纠正过来。
假设f(x,y)是一个离散空间中的二维函数,则该函数的二维傅立叶变换的定义如下:
p=0,1…M-1 q=0,1…N-1 (1)
或 p=0,1…M-1 q=0,1…N-1 (2)
离散傅立叶反变换的定义如下:
m=0,1…M-1 n=0,1…N-1(3)
F(p,q)称为f(m,n)的离散傅立叶变换系数。这个式子表明,函数f(m,n)可以用无数个不同频率的复指数信号和表示,而在频率(w1,w2)处的复指数信号的幅度和相位是F(w1,w2)。
2、MATLAB提供的快速傅立叶变换函数
(1)fft2
fft2函数用于计算二维快速傅立叶变换,其语法格式为:
B = fft2(I)
B = fft2(I)返回图象I的二维fft变换矩阵,输入图象I和输出图象B大小相同。
例如,计算图象的二维傅立叶变换,并显示其幅值的结果,如图所示,其命令格式如下
load imdemos saturn2
imshow(saturn2)
B = fftshift(fft2(saturn2));
imshow(log(abs(B)),[],'notruesize')
(2)fftshift
MATLAB提供的fftshift函数用于将变换后的图象频谱中心从矩阵的原点移到矩阵的中心,其语法格式为:
B = fftshift(I)
对于矩阵I,B = fftshift(I)将I的一、三象限和二、四象限进行互换。
(2)ifft2
ifft2函数用于计算图象的二维傅立叶反变换,其语法格式为:
B = ifft2(I)
B = ifft2(A)返回图象I的二维傅立叶反变换矩阵,输入图象I和输出图象B大小相同。其语法格式含义与fft2函数的语法格式相同,可以参考fft2函数的说明。
如果信号是二维的,用上面的函数即可!直接调用。
如果信号是一维的,给你下面的例子,你应该能明白!
clear
fs=100;N=128; %采样频率和数据点数
n=0:N-1;t=n/fs; %时间序列
x=05sin(2pi15t)+2sin(2pi40t); %信号
y=fft(x,N); %对信号进行快速Fourier变换,逆变换函数为ifft
mag=abs(y); %求得Fourier变换后的振幅
f=nfs/N; %频率序列
subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=128');grid on;
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=128');grid on;
%对信号采样数据为1024点的处理
fs=100;N=1024;n=0:N-1;t=n/fs;
x=05sin(2pi15t)+2sin(2pi40t); %信号
y=fft(x,N); %对信号进行快速Fourier变换
mag=abs(y); %求取Fourier变换的振幅
f=nfs/N;
subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');grid on;
subplot(2,2,4)
plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');grid on;
以上就是关于做8点的FFT程序 用C或matlab实现全部的内容,包括:做8点的FFT程序 用C或matlab实现、如何使用fft函数进行编程序和进行快速傅里叶逆变换、等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!
欢迎分享,转载请注明来源:内存溢出
微信扫一扫
支付宝扫一扫
评论列表(0条)