matlab - 通过傅立叶空间填充进行插值

Question

我最近尝试在 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 的评论。

提前致谢

Answer 1

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

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

折叠 for 循环也是正确实现 matlab 的必要步骤，我不将其用于培训目的以简化我的理解和边界检查。

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

实际上我认为它在字面上等同于带有别名 sinc 函数的卷积，也称为狄利克雷核，当采样频率增加到无穷大时，它限制了经典的 sinc 函数(参见 http://www.dsprelated.com/dspbooks/sasp/Rectangular_Window_I.html#sec:rectwinintro)

我没有在这里发布整个代码，但我的原始程序的目的是比较 dirichlet 核卷积公式，我在不同的框架中演示(使用傅里叶级数离散表达式的理论演示)，sinc 卷积 Whittaker–Shannon interpolation formula , 和零填充，所以我应该得到一个非常相似的结果。

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

如果您的信号没有带宽限制，或者关于采样率混叠，您将需要减少频谱混叠部分的贡献，这是通过在频域中使用频率滤波器 = 切趾窗口​​滤除它们来完成的，至极转化为特定的时域插值核。