通过傅立叶空间填充进行插值 [英] Interpolation through fourier space padding

问题描述

我最近尝试在matlab上实现在傅立叶域中使用zéro填充的插值方法的简单示例. 但是我无法正常工作，我总是有一个很小的频移，在傅立叶空间中几乎看不到，但是在时空中会产生巨大的误差.

I recently tried to implement on matlab a simple example of interpolation method using zéro padding in the fourier domain. But I am not able to get this work properly, I always have a small frequency shift, barely not visible in fourier space, but thay generates a huge error in time space.

由于傅立叶空间中的zéro填充似乎是一种常见的(且快速的)插值方法，因此我认为我缺少一些东西:

As zéro padding in the fourier space seems to be a common (and fast) interpolation method, I assume that there is something I am missing:

这是matlab代码:

Here is the matlab code:

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的评论，我将不胜感激.

I am not able to post images of the result because of my reputation, but, graph are easy to generates, through matlab. I would really appreciate a comment on either this code or zero padding fft for interpolation in general.

提前谢谢

推荐答案

非常感谢你们两个人的建议，我决定回答我自己的问题，因为注释框内的可用空间不足:

Thank you very much for both of you for those advices, I decided to respond my own question because the was not enough space available in coment box:

@尽力而为，我确实定义了错误的离散时间向量，就像@Frederick指出的那样，我在填充向量时遇到了问题，感谢您为我提供了正确的"matlab方法"，我不应该这样做如此怕fftshift/ifftshift，此外，将padarray与'both'一起使用也可以完成工作，如@Frederick所提到的.

@Try hard I indeed defined a wrong discrete Time vector, as @Frederick also pointed out, I had a problem in padding my vector, thank you for giving me the right "matlab way" to do it, I should not have been so afraid of fftshift/ifftshift, in addition, the use of padarray with 'both' would also have done the job, as mentionned by @Frederick .

折叠for循环也是正确实施matlab的必不可少的步骤，我没有在训练中使用它来简化理解和边界检查.

Collapsing the for loop was also an essential step for proper matlab implementation, that I don't use in training purpose to ease my understanding and bound checking.

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

An additionnal very interesting point @Try hard mentionned in its first sentence, and that I did not realised in the first place, is the fact, that zero padding is just the equivalent of convoluting my data in time domain by a sinc function.

实际上，我认为它在某种程度上等同于具有别名sinc函数的卷积，也称为dirichlet核，当采样频率向无穷大增加时，它是经典的sinc函数(请参阅

Actually I think that it is literraly equivalent to a convolution with an aliased sinc function, also called dirichlet kernel, which limits, when sampling frequency increase towards infinity, is the classic sinc function (see http://www.dsprelated.com/dspbooks/sasp/Rectangular_Window_I.html#sec:rectwinintro)

我没有在此处发布整个代码，但是我的原始程序的目的是比较dirichlet内核卷积公式，该公式是我在不同的框架中演示的(使用傅立叶级数离散表达式的理论演示)，sinc卷积 Whittaker-Shannon插值公式，以及零填充，所以应该给我一个非常相似的结果.

I did not posted the whole code here, but the purpose of my original program was to compare dirichlet kernel convolution formula, that I demonstrated in a different framework (theoretical demonstration using fourier series discrete expression) , sinc convolution Whittaker–Shannon interpolation formula, and zero padding, so I should be given with a very similar result.

对于变迹问题，我认为真正的答案是，如果您的信号是有限带宽的，那么除了矩形窗口之外，您不需要其他变迹功能.

To the apodization question, I think that the true answer is that, if your signal is bandlimited, you don't need other apodization function than rectangular window.

如果您的信号没有带宽限制，或者在采样率方面没有混叠，则需要减少频谱混叠部分的贡献，方法是使用频率滤波器将其滤出=频域中的变迹窗口，在时域中变成特定的插值内核.

If your signal is not bandlimited, or aliased regarding the sampling rate, you will need to reduce the contribution of aliased part of the spectrum, which is done by filtering them out with a frequency filter = apodizing window in frequency domain, wich turns into a specific interpolation kernels in time domain.

这篇关于通过傅立叶空间填充进行插值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！