傅里叶变换中的假频**
所有波都可由间谐波叠加而来,傅里叶变换实质上是对某一复杂波的组成成分进行分析,即将复杂波分离为不同频率的间谐波:
这里以一简单信号作为示例:
图一信号可由图二、图三与图四叠加而来,图二、三、四均为简单间谐波,其频率分别为:20Hz、60Hz、85Hz。傅里叶变换即为叠加过程的逆过程。
在对某一信号做傅里叶变换时,需要格外注意采样间距dt,采样频率fn=1/dt,以及奈奎斯特频率(折叠频率)fs=fn/2,为了防止假频的出现,采样定理规定采样频率需大于信号最大频率的2倍,即fn>2fmax;而当采样频率 fn<2fmax时,假频可由fa=|mfn-f|求得,其中m为整数,这相当于将超出奈奎斯特频率的部分折叠回[-fs,fs]区间,导致假频的出现。下面上述信号作为例子:
上图第一张图为采样间距dt=1ms时对应频谱,此时其折叠频率fs为1/2*dt,即500Hz,此时无假频出现;
第二张图为采样间距dt=10ms时对应频谱,此时其折叠频率fs为50Hz,可观察到60Hz信号此时对应频率60-100=-40Hz,85Hz信号此时对应频率85-100=-15Hz;
第三张图为采样间距dt=20ms时对应频谱,此时其折叠频率fs为25Hz,可观察到60Hz信号此时对应频率60-50=10Hz,85Hz信号此时对应频率85-2*50=-15Hz。
以下为Matlab代码:
dt1=0.001;
dt2=0.01;
dt3=0.02;
f1=20;
f2=60;
f3=85;
t_start=-0.5;
t_end=0.5;
df=1/(t_end-t_start); %频率间隔
t1=t_start:dt1:t_end;
t2=t_start:dt2:t_end;
t3=t_start:dt3:t_end;
f1_real=zeros(1,(t_end-t_start)/dt1);
f1_image=f1_real;
f1_a=f1_real;
f2_real=zeros(1,(t_end-t_start)/dt2);
f2_image=f2_real;
f2_a=f2_real;
f3_real=zeros(1,(t_end-t_start)/dt3);
f3_image=f3_real;
f3_a=f3_real;
source1=sin(2pif1t1+pi/6)+sin(2pif2t1+pi/4)+sin(2pif3t1+pi/5);
source2=sin(2pif1t2+pi/6)+sin(2pif2t2+pi/4)+sin(2pif3t2+pi/5);
source3=sin(2pif1t3+pi/6)+sin(2pif2t3+pi/4)+sin(2pif3t3+pi/5);
% plot(t,source);
for j=1:df:1/dt1
for i=t_start:dt1:t_end-dt1
z=floor((i-t_start)/dt1+1);
f1_real(j)=f1_real(j)+cos(2piji)source1(z);
f1_image(j)=f1_image(j)+sin(2piji)*source1(z);
end
f1_a(j)=sqrt(f1_real(j)*f1_real(j)+f1_image(j)*f1_image(j));
end
for j=1:df:1/dt2
for i=t_start:dt2:t_end-dt2
z=floor((i-t_start)/dt2+1);
f2_real(j)=f2_real(j)+cos(2piji)source2(z);
f2_image(j)=f2_image(j)+sin(2pij*i)*source2(z);
end
f2_a(j)=sqrt(f2_real(j)*f2_real(j)+f2_image(j)*f2_image(j));
end
for j=1:df:1/dt3
for i=t_start:dt3:t_end-dt3
z=floor((i-t_start)/dt3+1);
f3_real(j)=f3_real(j)+cos(2piji)source3(z);
f3_image(j)=f3_image(j)+sin(2pij*i)*source3(z);
end
f3_a(j)=sqrt(f3_real(j)*f3_real(j)+f3_image(j)*f3_image(j));
end
% f2=fft(source2);
figure(1)
subplot(3,1,1);
% plot(-1/dt1/2:df:1/dt1/2-1,fftshift(f1_a));
plot(-1/dt1/2:df:1/dt1/2,abs(fftshift(fft(source1))));
subplot(3,1,2);
plot(-1/dt2/2:df:1/dt2/2,abs(fftshift(fft(source2))));
% plot(-1/dt2/2:df:1/dt2/2-1,abs(fftshift(f2_a)));
subplot(3,1,3)
plot(-1/dt3/2:df:1/dt3/2,abs(fftshift(fft(source3))));
% plot(-1/dt3/2:df:1/dt3/2-1,abs(fftshift(f3_a)));
实际上该代码并没有利用sin与cos做加法的部分,我对比了求和与直接利用fft函数求出的结果,发现使用求和得出的结果与f使用ft函数得出结果基本一致,但频谱波折很大,暂时未找到原因,就先挂在这等以后再找一找问题。