3
我最近試圖在matlab上實現一個簡單的插值方法的例子，在傅立葉域中使用zéro填充。 但是我無法正確地完成這項工作，我總是有一個小的頻率偏移，在國際空間中幾乎看不到，但是它在時間空間中產生了一個巨大的錯誤。通過fourier空間填充插值
如傅立葉空間零填充似乎是一個常見（快速）插值方法，我認爲有我丟失的東西：
這裏是MATLAB代碼：
clc;
clear all;
close all;
Fe = 3250;
Te = 1/Fe;
Nech = 100;
F1 = 500;
F2 = 1000;
FMax = 1500;
time = [0:Te:(Nech-1)*Te];
timeDiscrete = [1:1:Nech];
frequency = (timeDiscrete/Nech)*Fe;
signal = cos(2*pi*F1*(time))+cos(2*pi*F2*(time))+cos(2*pi*FMax*(time));
%Compute the FFT
spectrum=zeros(1,Nech);
for k = timeDiscrete
for l = timeDiscrete
spectrum(k) = spectrum(k) + signal(l)*exp(-2*pi*j*l*k/Nech);
end
end
%Compute de inverse FFT
reconstruction=zeros(1,Nech);
for k = timeDiscrete
for l = timeDiscrete
reconstruction(k) = reconstruction(k) + spectrum(l)*exp(2*pi*j*l*k/Nech);
end
end
reconstruction=reconstruction/Nech;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% Now interpolation will take place %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Finterp = 6*Fe;
Tinterp = 1/Finterp;
TimeInterp = [0:Tinterp:(Nech-1)*Te];
[m,n] = size(TimeInterp);
NechInterp = n;
TimeInterpDiscrete = [1:1:NechInterp];
%Compute original signal value without any interpolation
signalResampled = cos(2*pi*F1*(TimeInterp))+cos(2*pi*F2*(TimeInterp))+cos(2*pi*FMax*(TimeInterp));
%Compute original signal interpolation by padding the fft and performing
%inverse fft on the result
semipaddedsize=floor(NechInterp/2);
padded_spectrum0 = zeros(1,semipaddedsize);
padded_spectrum0 = padarray(spectrum(1:Nech/2),[0 semipaddedsize-(Nech/2)],0,'post');
padded_spectrum = zeros(1,NechInterp);
padded_spectrum(1:semipaddedsize) = padded_spectrum0;
padded_spectrum(semipaddedsize+1:NechInterp-1) = conj(fliplr(padded_spectrum0));
% padded_spectrum = padarray(spectrum,[0 NechInterp-Nech],0,'post');
padded_timeDiscrete = [1:1:NechInterp];
padded_reconstruction = zeros(1,NechInterp);
for k = padded_timeDiscrete
for l = padded_timeDiscrete
padded_reconstruction(k) = padded_reconstruction(k) + padded_spectrum(l)*exp(2*pi*j*l*k/NechInterp);
end
end
padded_reconstruction=padded_reconstruction/(1*Nech);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% Let's print out the result %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
spectrumresampled=zeros(1,NechInterp);
for k = TimeInterpDiscrete
for l = TimeInterpDiscrete
spectrumresampled(k) = spectrumresampled(k) + signalResampled(l)*exp(-2*pi*j*l*k/NechInterp);
end
end
figure(2);
plot(abs(spectrumresampled)/6,'g');
hold on;
plot(abs(padded_spectrum),'b');
figure(3);
% Ground truth : deterministic signal is recomputed
plot(TimeInterp,signalResampled,'g');
hold on;
% linear interpolation between subsampled points (matlab tracing tool)
plot(time,(reconstruction),'c');
hold on;
% Padding the spectrum and reconstructing it
plot(TimeInterp,real(padded_reconstruction),'m');
hold on;
xlabel('Time in s','FontSize',16)
ylabel('Signal value (no unit)','FontSize',16)
title('\it{ Various signal reconstruction from fourier transform }','FontSize',16)
legend('True signal', 'Reconstruction with linear interpolation', 'Reconstruction with padded spectrum');
由於我的聲望，我無法發佈結果圖像，但通過matlab，圖形很容易生成。 我真的很感激對這個代碼或零填充fft一般插值的評論。
預先感謝您
你知道嗎？唯一的'fft'不能代表你信號的功率譜？ – fpe
其實我這裏不使用FFT，我用手寫的DFT。我知道我可以通過計算得到的dft完美地重建我的原始信號，當我修改頻譜以獲得更好的時間空間分辨率時，會出現問題。據我所知，許多應用程序正在使用fft和ifft來執行快速和準確的插值，並且我從來沒有聽說過會導致問題的傅立葉變換或其實現的理論侷限性？你指的是什麼？ – Tobbey
Tobbey，我認爲你的DFT有些問題。只運行第一塊代碼，然後執行：繪製信號;等一下;情節（ABS（重建）， 'G'）;顯示我他們不排隊，並且一個被翻轉和不同的幅度... – Frederick