clc
clear
close all,
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 200; % Length of signal
t = (0:L-1)*T; % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180);
y = x+2*randn(size(t)); % Sinusoids plus noise
figure,
NFFT = 2^nextpow2(L); % Next power of 2 from length of y比L大的
Y = fft(y,NFFT)/L; %要除以原始信号的采样个数数才能得到准确的振幅,
f = Fs/2*linspace(0,1,NFFT/2+1);%对于1~256个点,其对称点出现在129处
% Plot single-sided amplitude spectrum.
plot(f,2.*abs(Y(1:NFFT/2+1))) %且对于物理频谱即单边谱则需取其中一半,然后乘2
%分配到0—Fs之间,注意是将NFFT个数分配到此区间,即NFFT-1等分。...
%对称是以除去第一个点后剩下的奇数个数的中点为对称轴的,且对称是指共轭对称,0频率的振幅为0
title('单边谱')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
figure,
f = Fs*linspace(0,1,NFFT);
plot(abs(Y(1:NFFT)),'r')
title('双边谱')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
figure,
%双边谱恢复
yi=ifft(Y).*L;%完全恢复,针对傅氏变换的双边谱复数
yi=real(yi(1:L));
plot(t,yi)
hold on,
plot(t,y,'r')
title('双边谱恢复')
legend('恢复信号','真信号')
figure,
%单边谱恢复
for i=1:length(Y)
if(i<=length(Y)/2+1)
Y_d(i)=2*Y(i);%取一半后,能量守恒,振幅就要加倍
else
Y_d(i)=0;%将单边谱缺的那一段补0
end
end
yii=ifft(Y_d).*L;%完全恢复,要乘以L原始信号采样个数,因为之前除了采样个数
yii=real(yi(1:L));%对于fft其频谱为复数的模,ifft其恢复信号为复数的实部
plot(t,yii)
hold on,
plot(t,y,'r')
title('单边谱恢复')
legend('恢复信号','真信号')
%%%%对频率域的对称复数列的两边都要做相应的处理,再反傅氏变换才能回去,
%%%也可以只对1:129的做衰减,将130-256补零,其后的数在反傅氏变换后只影响虚部
%%%而实部才是我们要的,但此时需要将恢复的振幅乘以2,因为恢复时少了一半的能量
%%复信号代表相位不为零
%%频率域必须以复数的模的2倍作为真实振幅谱,反变换是必须是2^n个的复数
%%%%%结论:时间域数列补零到2^n,对频率域的复数列也为2^n,...
%直接对其做直求反变换再截断到原始时间域的长度即为最终的恢复数列(若返回的是复数,则取其实部即可)
%tips:使用基2算法时,必须将数补零为2的n次方,此时参与fft运算的数据为2^n个,
%得到的结果也为2^n个,然后要将能量均分,即除以2^n,保证能量守恒