通过傅里叶空间填充进行插值

我最近尝试在 matlab 上实现一个在傅立叶域中使用零填充的插值方法的简单示例。 但我无法正常工作，我总是有一个小的频移，在傅里叶空间中几乎不可见，但它在时空上产生了巨大的误差。

由于傅里叶空间中的零填充似乎是一种常见（且快速）的插值方法，因此我认为我缺少一些东西：

这是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 的评论。

先感谢您

非常感谢你们两位的建议，我决定回答我自己的问题，因为评论框中没有足够的空间：

@努力尝试我确实定义了一个错误的离散时间向量，正如@Frederick也指出的那样，我在填充我的向量时遇到了问题，谢谢你给我正确的“matlab方式”来做到这一点，我不应该这么害怕fftshift/ifftshift，此外，使用 padarray 和“both”也可以完成这项工作，正如 @Frederick 提到的。

折叠 for 循环也是正确实施 matlab 的重要步骤，我不会在训练中使用它来简化我的理解和边界检查。

@Try Hard 在第一句话中提到的另一个非常有趣的点是我一开始没有意识到的事实是，零填充相当于通过 sinc 函数在时域中对我的数据进行卷积。

实际上，我认为它在字面上相当于具有别名 sinc 函数的卷积，也称为 dirichlet 内核，当采样频率增加到无穷大时，它限制了经典的 sinc 函数（请参阅http://www.dsplated.com/dspbooks/sasp/Rectangle_Window_I.html#sec:rectwinintro http://www.dsprelated.com/dspbooks/sasp/Rectangular_Window_I.html#sec%3arectwinintro)

我没有在这里发布完整的代码，但我原始程序的目的是比较狄利克雷核卷积公式，我在不同的框架中演示了该公式（使用傅立叶级数离散表达式的理论演示），sinc卷积Whittaker–Shannon 插值公式 http://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula和零填充，所以我应该得到非常相似的结果。

对于切迹问题，我认为真正的答案是，如果您的信号是带限的，则除了矩形窗口之外，您不需要其他切迹函数。

如果您的信号不受带宽限制，或者在采样率方面存在混叠，则需要减少频谱混叠部分的贡献，这是通过使用频率滤波器 = 频域中的变迹窗口将其过滤掉来完成的，这将变成时域中的特定插值内核。