Интерполация чрез запълване на пространството на Фурие

Благодаря ви много и на двама ви за тези съвети, реших да отговоря на собствения си въпрос, защото нямаше достатъчно място в полето за коментари:

@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)

Не съм публикувал целия код тук, но целта на оригиналната ми програма беше да сравня формулата за навиване на ядрото на Дирихле, която демонстрирах в различна рамка (теоретична демонстрация, използваща дискретни изрази от серия на Фурие), sinc навивка Интерполационна формула на Whittaker–Shannon и нулево подпълване, така че трябва да бъдат дадени с много подобен резултат.

На въпроса за аподизация, мисля, че истинският отговор е, че ако сигналът ви е с ограничен обхват, нямате нужда от друга функция за аподизация освен правоъгълен прозорец.

Ако вашият сигнал не е с ограничена честотна лента или псевдоним по отношение на честотата на дискретизация, ще трябва да намалите приноса на част от спектъра с псевдоним, което се прави чрез филтрирането им с честотен филтър = аподизиращ прозорец в честотната област, който се превръща в специфични интерполационни ядра във времева област.

ДОБРЕ. Един проблем беше начинът, по който правите IDFT за padded_reconstruction. Начинът, по който дефинирахте TimeInterp и по този начин NechInterp направиха елементите на сложния показател неправилни. Това обяснява неправилните честоти.

След това имаше проблем с включването на средната точка в домейна на Фурие (т. 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);

Резултатът, който наблюдавате във времевата област, е звънене поради конволюция на sinc функция с вашите оригинални данни. Това е еквивалентът във времевата област на умножение с правоъгълен прозорец в честотната област, което всъщност е това, което правите, когато запълвате нула. Не забравяйте да apodize!

Публикувам повторно вашия код след свиване на циклите (което значително ускорява изчислението), предефиниране на обхватите на променливите за време и честота (вижте дефиницията на 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;

След това стигате до този много по-добър интерполиран сигнал във времева област, тъй като операцията за подплънка не води до такива резки спадове в спектъра (все пак може да са възможни някои подобрения с по-нататъшно бърникане):