3

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.

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:

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');

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 in advance

4

3 回答 3

2

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

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

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

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

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

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

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

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

于 2013-08-21T09:59:34.867 回答
1

由于 sinc 函数与原始数据的卷积，您在时域中观察到的结果正在振铃。这相当于在时域中乘以频域中的矩形窗口，这实际上是您在零填充时所做的事情。不要忘记变迹！

在折叠循环（显着加速计算）、重新定义时间和频率变量的范围（请参阅DFT的定义以了解原因）并删除您执行的其中一个填充操作后，我重新发布您的代码，坦率地说，我没有明白点了。

clc;
clear all;
close all;

Fe = 3250;
Te = 1/Fe;
Nech = 100;
mlt = 10;
F1 = 50; 
F2 = 100; 
FMax = 150; 

time = [0:Te:(Nech-1)*Te];
%timeDiscrete = [1:1:Nech];
timeDiscrete = [0:1:Nech-1];
frequency = (timeDiscrete/Nech)*Fe;

signal = cos(2*pi*F1*(time))+cos(2*pi*F2*(time))+cos(2*pi*FMax*(time));
spectrum =  signal*exp(-2*pi*j*timeDiscrete'*timeDiscrete/Nech);
fspec = [0:Nech-1]*Fe/Nech;
reconstruction = spectrum*exp(2*pi*j*timeDiscrete'*timeDiscrete/Nech)/Nech;

figure
plot(time,signal)
hold on
plot(time,reconstruction,'g:')

% **** interpolation ****

Finterp = 6*Fe;
Tinterp = 1/Finterp;
TimeInterp = [0:Tinterp:(Nech-1)*Te];
NechInterp = length(TimeInterp);
%TimeInterpDiscrete = [1:NechInterp];
TimeInterpDiscrete = [0:NechInterp-1];

%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
padded_spectrum0 = spectrum;
padded_spectrum0(NechInterp) = 0;
fspecPadded = [0:NechInterp-1]*Finterp/NechInterp;

padded_reconstruction = padded_spectrum0*exp(2*pi*j*TimeInterpDiscrete'*TimeInterpDiscrete/NechInterp)/(1*Nech);
spectrumresampled = signalResampled*exp(-2*pi*j*TimeInterpDiscrete'*TimeInterpDiscrete/NechInterp);
fresampled = [0:NechInterp-1]*Fe/NechInterp;

% **** print out ****

figure(2);
hold on;
plot(fspec,abs(spectrum),'c');
plot(fresampled,abs(spectrumresampled)/6,'g--');
plot(fspecPadded,abs(padded_spectrum0),'m:');
xlabel('Frequency in Hz','FontSize',16)
ylabel('Signal value (no unit)','FontSize',16)
legend('True signal', 'Reconstruction with resampled spectrum', 'Padded spectrum');

figure(3);
% Ground truth : deterministic signal is recomputed
plot(TimeInterp,signalResampled,'g');
hold on;
% Padding the spectrum and reconstructing it 
plot(TimeInterp,real(padded_reconstruction),'m:');
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 padded spectrum');

这是由于频域中的填充而导致的严重失真信号的图像：

在此处输入图像描述

通过首先将fftshift频谱居中并在居中频谱的另一侧进行填充，然后反转fftshift操作，可以进行一些改进：

Nz = NechInterp-Nech;
padded_spectrum0 = ifftshift([ zeros(1,floor(Nz/2)) fftshift(spectrum) zeros(1,floor(Nz/2)+rem(Nz,2)) ]); % replaces (NechInterp) = 0;
fspecPadded = [0:NechInterp-1]*Finterp/NechInterp;

然后你会得到这个更好的插值时域信号，因为填充操作不会导致频谱出现如此突然的下降（一些改进可能仍然是可能的，但需要进一步修补）：

在此处输入图像描述

于 2013-08-20T22:42:34.763 回答
1

好的。一个问题是您为 padded_reconstruction 执行 IDFT 的方式。您定义 TimeInterp 并因此定义 NechInterp 的方式使复数指数的元素不正确。这解释了不正确的频率。

然后在傅立叶域（pt 50）中包含两次中点存在问题。它接近于零，所以它不是一个非常明显的问题，但它应该只包含一次。我重写了这一部分，因为我很难完全按照你的方式完成它。不过，我把它保持得很近。如果我这样做，我会使用 fftshift，然后使用 padarray(...,'both')，这样可以省去将零放在中间的工作。如果您这样做是为了获得学习经验并尝试不使用 matlab 工具（例如 fftshift），那么没关系。

我还重新定义了你定义时间的方式，但公平地说，我认为它可以按照你的方式工作。

我已经用 %<<<<<<<<<<

Fe = 3250;
Te = 1/Fe;
Nech = 100;

F1 = 500;
F2 = 1000;
FMax = 1500;

time = [Te:Te:(Nech)*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 = [Tinterp:Tinterp:(Nech*6)*Tinterp];  %<<<<<<<<<<
[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
padded_spectrum = zeros(1,NechInterp); %<<<<<<<<<<
padded_spectrum(1:floor(Nech/2-1)) = spectrum(1:floor(Nech/2-1)); %<<<<<<<<<<
padded_spectrum(end-floor(Nech/2)+1:end) = spectrum(floor(Nech/2)+1:end); %<<<<<<<<<<

padded_reconstruction = zeros(1,NechInterp);
for k = TimeInterpDiscrete %<<<<<<<<<<(no reason for new variable)
    for l = TimeInterpDiscrete %<<<<<<<<<<(no reason for new variable)
        padded_reconstruction(k) = padded_reconstruction(k) + padded_spectrum(l)*exp(2*pi*j*l*k/NechInterp);
    end
end
padded_reconstruction=padded_reconstruction/(1*Nech);
于 2013-08-20T22:44:38.783 回答