matlab中如何对一组复数频域信号进行快速傅里叶逆变换
matlab中如何对一组复数频域信号进行快速傅里叶逆变换为何需要对复数频域信号进行快速傅里叶逆变换快速傅里叶变换(fft)后的数据格式快速傅里叶逆变换至原信号时域波型第一步第二步第三步为何需要对复数频域信号进行快速傅里叶逆变换目前查到很多资料对快速傅里叶变换讲解非常清楚,包括其原理以及变换后的具体处理均有很多介绍,可参见以下链接:但是对如何处理一组频域信号转换成原始信号的傅里叶逆变换很少介绍,并且
matlab中如何对一组复数频域信号进行快速傅里叶逆变换
为何需要对复数频域信号进行快速傅里叶逆变换
目前查到很多资料对快速傅里叶变换讲解非常清楚,包括其原理以及变换后的具体处理均有很多介绍,可参见以下链接:
但是对如何处理一组频域信号转换成原始信号的傅里叶逆变换很少介绍,并且直接对一组频域信号进行快速傅里叶逆变换(ifft)是还原不到原始时域状态的,经过一番研究后,在此进行相应的整理,以便日后进行使用。
快速傅里叶变换(fft)后的数据格式
matlab中一个信号经过快速傅里叶变换后,得到的数据与原始数据的采样个数是一致的,但是这里面有一半的数据是重复的,举例:
采样N=256个数据,采样频率是Fs=256,信号为
S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)
对其进行快速傅里叶变换
代码如下
clear;
Fs=256;
T=1/Fs;
N=256;
t=(0:N-1)*T;
S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180);
Y=fft(S);
Y即为快速傅里叶变换后的数据,这组数据的格式如下:
- 每个数据的实部为an,虚部位bn,在matlab中看到的形式是an+bni,sqrt(an2+bn2)为其模,arctan(b/a)为其相位角。
- Y(1)为直流分量的N倍,即0Hz时的值再乘以N,此例中Y(1)为512(即2*256);
- Y(2:N/2+1)为N/2倍的信号值,对应频率从Fs/N ~ Fs/2;
- Y(N/2+1:end)对应频率从-Fs/2 ~ -Fs/N,可以发现Y(N/2+1)点被共用,两值是相加,所以该点比较特殊。这部分的值和Y(2:N/2+1)的值是共轭的,即实部相同虚部相反,这也是为什么说快速傅里叶变换出来的值有一半是重复的原因。
此部分可参考链接: 关于实信号的双边谱和单边谱.
上图为对Y直接求模后的图形,横坐标为其对应的频率,可看到数据关于第N/2+1点(即第127点)对称,Y(1)点对应的频率为0Hz,值为2x256=512。
代码复现上图如下:
clear;
Fs=256;
T=1/Fs;
N=256;
Fn=zeros(N,1);
t=(0:N-1)*T;
S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180);
for n=1:N
Fn(n)=(n-1)*Fs/N;
end
Y=fft(S);
plot(Fn,abs(Y));
从以上数据结构中就不难发现,如果你直接对一个信号进行逆向变换时,它并不满足上面的这种格式,所以变化出来的值必定不是原信号的模样。所以为了能够复现出原来的信号,必须对现有信号按照以上格式进行修改才行。
快速傅里叶逆变换至原信号时域波型
聪明的你看到这应该知道怎么去做了,以下仅是举例说明。
假设这时你有一个频率信号D(x),x为频率,它的频谱范围为0-128Hz,以1Hz为频率间隔,总共129个数据点,并且每个点的数据为an+bni的形式存在。
第一步
保留直流分量到数据首部,即Y(1)=D(0)x128x2。
第二步
Y(2:128)=D(1:127)x128
Y(256:130)=D(1:127)x128
Y(129)=Real(D(128)x256)
第三步
对Y进行快速逆变换就可以得到原始的数据图形
细心的你会发现Y(128)中只取了实部进行计算,这是因为此点为正频与负频共用点,两个值是共轭的,所以他们相加后得到的值仅有实值,虚值部分抵消了。
验证
现有一组数据点击此处下载数据,如下图所示
该数据为[100,5]的参数,第一列为频率,第二列为幅值,第三列为相位角(角度单位),第四列为实值,第五列为虚值。
采用上文中的步骤对数据进行处理,代码如下:
clear;
data=load('test.txt');
[m,n]=size(data);
N=m*2;
predata=zeros(2*m,1);
predata(1)=0; %由于本数据未给出直流分量,所以默认为0
for i=2:N/2
predata(i)=(data(i-1,4)+data(i-1,5)*sqrt(-1))*(N/2);%对2~N/2的数据进行赋值
end
predata(N/2+1)=data(N/2,4)*N;%对第N/2+1位置进行赋值
for i=2:N/2
predata(N-i+2)=(data(i-1,4)-data(i-1,5)*sqrt(-1))*(N/2);%对N/2+2~END的数据进行赋值
end
%此时predata的数据已处理完成,直接做逆变换即可
y=real(ifft(predata));
%接下来需要对时间参量进行确定
%数据的采集点总数N已知,需要确认采样频率Fs,可求出采样周期T
Fn=data(2,1)-data(1,1);
Fs=Fn*N;
T=1/Fs;
for i=1:N
t(i)=(i-1)*T;
end
figure(1);
plot(t,y(:,1));
%{
根据傅里叶逆变换的定义,可以知道时域信号是由各种频域信号叠加而成,
那么我们根据这个特性直接将频域信号进行叠加来验证逆变换的正确性
%}
C(:)=data(:,2); %幅值
%相位角
phi=zeros(m);
nw=zeros(m);
for i=1:m
phi(i)=data(i,3)*pi/180;
nw(i)=data(i,1)*pi*2;
end
fz=zeros(m,N);
z=zeros(N,1);
for i=1:m
for j=1:N
fz(i,j)=C(i)*cos(nw(i)*t(j)+phi(i));
end
end
for i=1:N
z(i,1)=sum(fz(:,i));
end
figure(2);
plot(t,z);
figure(1) 经处理后的快速傅里叶逆变换图形
figure(2)利用傅里叶变换的数学意义进行变换后的图形
上两张图进行对比发现,波形和相位以及幅值基本一致(仅在负极大值处有偏差),说明数据经傅里叶还原后与原数据图像基本相同,还原成功!
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)