y=[18 2 22 32 35 37]
y1=fft(y)
结果:
y =
18000 20000 22000 32000 35000 37000
y1 =164000 -14000 + 25981i -07000 + 03464i -14000 -07000 - 03464i -14000 - 25981i
1、采样数据导入Matlab 。
采样数据的导入至少有三种方法。
第一就是手动将数据整理成Matlab支持的格式,这种方法仅适用于数据量比较小的采样。
第二种方法是使用Matlab的可视化交互 *** 作,具体 *** 作步骤为:File --> Import Data,然后在d出的对话框中找到保存采样数据的文件,根据提示一步一步即可将数据导入。这种方法适合于数据量较大,但又不是太大的数据。
第三种方法,使用文件读入命令。数据文件读入命令有textread、fscanf、load等,如采样数据保存在txt文件中,则推荐使用 textread命令。如[a,b]=textread('datatxt','%f%f%f'); 这条命令将datatxt中保存的数据三个三个分组,将每组的第一个数据送给列向量a,第三个数送给列向量b,第二个数据丢弃。命令类似于C语言,详细可查看其帮助文件。文件读入命令录入采样数据可以处理任意大小的数据量,且录入速度相当快,一百多万的数据不到20秒即可录入。
2、对采样数据进行频谱分析 。
频谱分析自然要使用快速傅里叶变换FFT了,对应的命令即 fft ,简单使用方法为:Y=fft(b,N),其中b即是采样数据,N为fft数据采样个数。一般不指定N,即简化为Y=fft(b)。Y即为FFT变换后得到的结果,与b的元素数相等,为复数。以频率为横坐标,Y数组每个元素的幅值为纵坐标,画图即得数据b的幅频特性;以频率为横坐标,Y数组每个元素的角度为纵坐标,画图即得数据b的相频特性。典型频谱分析M程序举例如下: clc fs=100;
t=[0:1/fs:100];
N=length(t)-1;%减1使N为偶数 %频率分辨率F=1/t=fs/N
p=13sin(0482pit)+21sin(0522pit)+11sin(0532pit) +05sin(182pit)+09sin(222pit);
%上面模拟对信号进行采样,得到采样数据p,下面对p进行频谱分析
figure(1) plot(t,p); grid on
title('信号 p(t)'); xlabel('t') ylabel('p') Y=fft(p);
magY=abs(Y(1:1:N/2))2/N; f=(0:N/2-1)'fs/N; figure(2)
%plot(f,magY);
h=stem(f,magY,'fill','--');
set(h,'MarkerEdgeColor','red','Marker','') grid on
title('频谱图 (理想值:[048Hz,13]、[052Hz,21]、[053Hz,11]、[18Hz,05]、[22Hz,09]) '); xlabel('f (Hz)') ylabel('幅值')
对于现实中的情况,采样频率fs一般都是由采样仪器决定的,即fs为一个给定的常数;另一方面,为了获得一定精度的频谱,对频率分辨率F有一个人为的规定,一般要求F<001,即采样时间ts>100秒;由采样时间ts和采样频率fs即可决定采样数据量,即采样总点数N=fsts。这就从理论上对采样时间ts和采样总点数N提出了要求,以保证频谱分析的精准度。
#include <mathh>
#include <stdioh>
#define N 8
void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il);
void main()
{
double xr[N],xi[N],Yr[N],Yi[N],l=0,il=0;
int i,j,n=N,k=3;
for(i=0;i<N;i++)
{
xr[i]=i;
xi[i]=0;
}
printf("------FFT------\n");
l=0;
kkfft(xr,xi,n,k,Yr,Yi,l,il);
for(i=0;i<N;i++)
{
printf("%-11lf + j %-11lf\n",Yr[i],Yi[i]);
}
printf("-----DFFT-------\n");
l=1;
kkfft(Yr,Yi,n,k,xr,xi,l,il);
for(i=0;i<N;i++)
{
printf("%-11lf + j %-11lf\n",xr[i],xi[i]);
}
getch();
}
void kkfft(double pr[], double pi[], int n, int k, double fr[], double fi[], int l, int il)
{
int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;
for (it=0; it<=n-1; it++)
{
m = it;
is = 0;
for(i=0; i<=k-1; i++)
{
j = m/2;
is = 2is+(m-2j);
m = j;
}
fr[it] = pr[is];
fi[it] = pi[is];
}
pr[0] = 10;
pi[0] = 00;
p = 6283185306/(10n);
pr[1] = cos(p);
pi[1] = -sin(p);
if (l!=0)
pi[1]=-pi[1];
for (i=2; i<=n-1; i++)
{
p = pr[i-1]pr[1];
q = pi[i-1]pi[1];
s = (pr[i-1]+pi[i-1])(pr[1]+pi[1]);
pr[i] = p-q;
pi[i] = s-p-q;
}
for (it=0; it<=n-2; it=it+2)
{
vr = fr[it];
vi = fi[it];
fr[it] = vr+fr[it+1];
fi[it] = vi+fi[it+1];
fr[it+1] = vr-fr[it+1];
fi[it+1] = vi-fi[it+1];
}
m = n/2;
nv = 2;
for (l0=k-2; l0>=0; l0--)
{
m = m/2;
nv = 2nv;
for(it=0; it<=(m-1)nv; it=it+nv)
for (j=0; j<=(nv/2)-1; j++)
{
p = pr[mj]fr[it+j+nv/2];
q = pi[mj]fi[it+j+nv/2];
s = pr[mj]+pi[mj];
s = s(fr[it+j+nv/2]+fi[it+j+nv/2]);
poddr = p-q;
poddi = s-p-q;
fr[it+j+nv/2] = fr[it+j]-poddr;
fi[it+j+nv/2] = fi[it+j]-poddi;
fr[it+j] = fr[it+j]+poddr;
fi[it+j] = fi[it+j]+poddi;
}
}
/逆傅立叶变换/
if(l!=0)
{
for(i=0; i<=n-1; i++)
{
fr[i] = fr[i]/(10n);
fi[i] = fi[i]/(10n);
}
}
/是否计算模和相角/
if(il!=0)
{
for(i=0; i<=n-1; i++)
{
pr[i] = sqrt(fr[i]fr[i]+fi[i]fi[i]);
if(fabs(fr[i])<0000001fabs(fi[i]))
{
if ((fi[i]fr[i])>0)
pi[i] = 900;
else
pi[i] = -900;
}
else
pi[i] = atan(fi[i]/fr[i])3600/6283185306;
}
}
return;
}
快速傅里叶变换 要用C++ 才行吧 你可以用MATLAB来实现更方便点啊
此FFT 是用VC60编写,由FFTCPP;STDAFXH和STDAFXCPP三个文件组成,编译成功。程序可以用文件输入和输出为文件。文件格式为TXT文件。测试结果如下:
输入文件:8TXT 或手动输入
8 //N
1
2
3
4
5
6
7
8
输出结果为:或保存为TXT文件。(8OUTTXT)
8
(36,0)
(-4,965685)
(-4,4)
(-4,165685)
(-4,0)
(-4,-165685)
(-4,-4)
(-4,-965685)
下面为FFTCPP文件:
// FFTcpp : 定义控制台应用程序的入口点。
#include "stdafxh"
#include <iostream>
#include <complex>
#include <bitset>
#include <vector>
#include <conioh>
#include <string>
#include <fstream>
using namespace std;
bool inputData(unsigned long &, vector<complex<double> >&); //手工输入数据
void FFT(unsigned long &, vector<complex<double> >&); //FFT变换
void display(unsigned long &, vector<complex<double> >&); //显示结果
bool readDataFromFile(unsigned long &, vector<complex<double> >&); //从文件中读取数据
bool saveResultToFile(unsigned long &, vector<complex<double> >&); //保存结果至文件中
const double PI = 31415926;
int _tmain(int argc, _TCHAR argv[])
{
vector<complex<double> > vecList; //有限长序列
unsigned long ulN = 0; //N
char chChoose = ' '; //功能选择
//功能循环
while(chChoose != 'Q' && chChoose != 'q')
{
//显示选择项
cout << "\nPlease chose a function" << endl;
cout << "\t1Input data manually, press 'M':" << endl;
cout << "\t2Read data from file, press 'F':" << endl;
cout << "\t3Quit, press 'Q'" << endl;
cout << "Please chose:";
//输入选择
chChoose = getch();
//判断
switch(chChoose)
{
case 'm': //手工输入数据
case 'M':
if(inputData(ulN, vecList))
{
FFT(ulN, vecList);
display(ulN, vecList);
saveResultToFile(ulN, vecList);
}
break;
case 'f': //从文档读取数据
case 'F':
if(readDataFromFile(ulN, vecList))
{
FFT(ulN, vecList);
display(ulN, vecList);
saveResultToFile(ulN, vecList);
}
break;
}
}
return 0;
}
bool Is2Power(unsigned long ul) //判断是否是2的整数次幂
{
if(ul < 2)
return false;
while( ul > 1 )
{
if( ul % 2 )
return false;
ul /= 2;
}
return true;
}
bool inputData(unsigned long & ulN, vector<complex<double> >& vecList)
{
//题目
cout<< "\n\n\n==============================Input Data===============================" << endl;
//输入N
cout<< "\nInput N:";
cin>>ulN;
if(!Is2Power(ulN)) //验证N的有效性
{
cout<< "N is invalid (N must like 2, 4, 8, ), please retry" << endl;
return false;
}
//输入各元素
vecListclear(); //清空原有序列
complex<double> c;
for(unsigned long i = 0; i < ulN; i++)
{
cout << "Input x(" << i << "):";
cin >> c;
vecListpush_back(c);
}
return true;
}
bool readDataFromFile(unsigned long & ulN, vector<complex<double> >& vecList) //从文件中读取数据
{
//题目
cout<< "\n\n\n===============Read Data From File==============" << endl;
//输入文件名
string strfilename;
cout << "Input filename:" ;
cin >> strfilename;
//打开文件
cout << "open file " << strfilename << "" <<endl;
ifstream loadfile;
loadfileopen(strfilenamec_str());
if(!loadfile)
{
cout << "\tfailed" << endl;
return false;
}
else
{
cout << "\tsucceed" << endl;
}
vecListclear();
//读取N
loadfile >> ulN;
if(!loadfile)
{
cout << "can't get N" << endl;
return false;
}
else
{
cout << "N = " << ulN << endl;
}
//读取元素
complex<double> c;
for(unsigned long i = 0; i < ulN; i++)
{
loadfile >> c;
if(!loadfile)
{
cout << "can't get enough infomation" << endl;
return false;
}
else
cout << "x(" << i << ") = " << c << endl;
vecListpush_back(c);
}
//关闭文件
loadfileclose();
return true;
}
bool saveResultToFile(unsigned long & ulN, vector<complex<double> >& vecList) //保存结果至文件中
{
//询问是否需要将结果保存至文件
char chChoose = ' ';
cout << "Do you want to save the result to file (y/n):";
chChoose = _getch();
if(chChoose != 'y' && chChoose != 'Y')
{
return true;
}
//输入文件名
string strfilename;
cout << "\nInput file name:" ;
cin >> strfilename;
cout << "Save result to file " << strfilename << "" << endl;
//打开文件
ofstream savefile(strfilenamec_str());
if(!savefile)
{
cout << "can't open file" << endl;
return false;
}
//写入N
savefile << ulN << endl;
//写入元素
for(vector<complex<double> >::iterator i = vecListbegin(); i < vecListend(); i++)
{
savefile << i << endl;
}
//写入完毕
cout << "save succeed" << endl;
//关闭文件
savefileclose();
return true;
}
void FFT(unsigned long & ulN, vector<complex<double> >& vecList)
{
//得到幂数
unsigned long ulPower = 0; //幂数
unsigned long ulN1 = ulN - 1;
while(ulN1 > 0)
{
ulPower++;
ulN1 /= 2;
}
//反序
bitset<sizeof(unsigned long) 8> bsIndex; //二进制容器
unsigned long ulIndex; //反转后的序号
unsigned long ulK;
for(unsigned long p = 0; p < ulN; p++)
{
ulIndex = 0;
ulK = 1;
bsIndex = bitset<sizeof(unsigned long) 8>(p);
for(unsigned long j = 0; j < ulPower; j++)
{
ulIndex += bsIndextest(ulPower - j - 1) ulK : 0;
ulK = 2;
}
if(ulIndex > p)
{
complex<double> c = vecList[p];
vecList[p] = vecList[ulIndex];
vecList[ulIndex] = c;
}
}
//计算旋转因子
vector<complex<double> > vecW;
for(unsigned long i = 0; i < ulN / 2; i++)
{
vecWpush_back(complex<double>(cos(2 i PI / ulN) , -1 sin(2 i PI / ulN)));
}
for(unsigned long m = 0; m < ulN / 2; m++)
{
cout<< "\nvW[" << m << "]=" << vecW[m];
}
//计算FFT
unsigned long ulGroupLength = 1; //段的长度
unsigned long ulHalfLength = 0; //段长度的一半
unsigned long ulGroupCount = 0; //段的数量
complex<double> cw; //WH(x)
complex<double> c1; //G(x) + WH(x)
complex<double> c2; //G(x) - WH(x)
for(unsigned long b = 0; b < ulPower; b++)
{
ulHalfLength = ulGroupLength;
ulGroupLength = 2;
for(unsigned long j = 0; j < ulN; j += ulGroupLength)
{
for(unsigned long k = 0; k < ulHalfLength; k++)
{
cw = vecW[k ulN / ulGroupLength] vecList[j + k + ulHalfLength];
c1 = vecList[j + k] + cw;
c2 = vecList[j + k] - cw;
vecList[j + k] = c1;
vecList[j + k + ulHalfLength] = c2;
}
}
}
}
void display(unsigned long & ulN, vector<complex<double> >& vecList)
{
cout << "\n\n===========================Display The Result=========================" << endl;
for(unsigned long d = 0; d < ulN;d++)
{
cout << "X(" << d << ")\t\t\t = " << vecList[d] << endl;
}
}
下面为STDAFXH文件:
// stdafxh : 标准系统包含文件的包含文件,
// 或是常用但不常更改的项目特定的包含文件
#pragma once
#include <iostream>
#include <tcharh>
// TODO: 在此处引用程序要求的附加头文件
下面为STDAFXCPP文件:
// stdafxcpp : 只包括标准包含文件的源文件
// FFTpch 将成为预编译头
// stdafxobj 将包含预编译类型信息
#include "stdafxh"
// TODO: 在 STDAFXH 中
//引用任何所需的附加头文件,而不是在此文件中引用
先给你个利用matlab中傅里叶变换进行函数频谱分析的程序。
clf;
fs=100;N=128; %采样频率和数据点数
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,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;
系统实现的意义和必要性:通过变换可以将原本的时域信号函数转化为频域,可以直观的观察到所采集信号函数的频域特征,更有利于进行信号分析
系统功能:通过傅里叶变换将信号函数的时域图形转化成频域图形,即将信号函数原本幅值随时间变化的特性曲线转化为幅值随频率变化的特性曲线
系统设计:利用傅里叶变换的快速傅里叶变换特性(fft)
系统测试:可以得出信号函数的功率谱函数,从时、频两域分析信号
遇到的问题及解决方法:mag=abs(y); 由于傅里叶变换是复数域的变换,需要对其函数进行求模
结论:利用傅里叶变换及其逆变换可以简单的将信号函数进行时、频域的转化,有利于进行信号分析。
本人有一定从事信号分析的经历,并且积累了一定的经验,对傅里叶变换有比较深的了解,希望对你有所帮助,如果有什么问题可以在线问我。
求xa=exp(-1000abs(t))在t=[-0005,0005]的傅里叶变换。
Dt=000005;
t=-0005:Dt:0005;
xa=exp(-1000abs(t)); %模拟信号
Wmax=2pi2000; %Dt=000005 so 周期为2pi2000
K=500;k=0:1:K;
W=kWmax/K; %将Wmax分为等间隔的500点,W是离散化后的旋转因子
Xa=xaexp(-jt'W)Dt;
Xa=real(Xa); %Xa=real(Xa)其实是取Xa各元素的模(幅值)
%连续时间傅立叶变换
W=[-fliplr(W),W(2:501)];
%频率从 -Wmax to Wmax
Xa=[fliplr(Xa), Xa(2:501)];
% Xa 范围 -Wmax to Wmax
figure(1)
subplot(2,1,1);
plot(t1000,xa,'');
xlabel('t in msec');
ylabel('xa(t)');
gtext('模拟信号');
subplot(2,1,2);
plot(W/(2pi1000),Xa1000,'');
xlabel('Frequence in KHz');
ylabel('Xa(jw)1000');
gtext('连续时间傅里叶变换');
以上就是关于在matlab中怎样用快速傅里叶变换求相位图 例如y=[1.8,2,2.2,3.2,3.5,3.7] 求程序全部的内容,包括:在matlab中怎样用快速傅里叶变换求相位图 例如y=[1.8,2,2.2,3.2,3.5,3.7] 求程序、matlab里有什么工具箱,可以用FFT(快速傅立叶变换)做频谱分析、傅里叶变换用C语言程序怎么实现等相关内容解答,如果想了解更多相关内容,可以关注我们,你们的支持是我们更新的动力!
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)