Idealna rekonstrukcja sygnału

Mamy sygnał fin, który został spróbkowany z częstotliwością fsin. Chcemy wygenerować sygnał, który został spróbkowany z częstotliwością fsout. Jeśli fsout > fsin w sygnale wyjściowym nie powinno być żadnej częstotliwości powyżej fsin. Jeśli fsout < fsin w sygnale wyjściowym powinny się znaleźć tylko częstotliwości mniejsze od fsout z sygnału wejściowego.

Wszelkie interpolacje (spline, bicubic, itp) w dziedzinie czasu odpadają, gdyż wprowadzimy tym częstotliwości do sygnału, których w nim nie było oryginalnie. Jedne sprawują się gorzej, inne lepiej, ale generalnie interpolacje odtworzą sygnał wizualnie niby dobrze, częstotliwościowo nie będzie to to samo. Ma to szczególne znaczenie kiedy pracujemy z sygnałem dźwiękowym. No dobra za wyjątkiem podanej niżej interpolacji. Rekonstrukcję możemy też przeprowadzić z dziedzinie częstotliwości.

Aby odtworzyć oryginalny sygnał potrzeba nam minimum dwie próbki na na okres, w praktyce powinno być ich więcej. Możemy próbkować dokładnie w miejscach gdzie sinusoida przechodzi przez zero. A nawet jeśli jesteśmy minimalnie minięci podczas próbkowania z zerem, pobierane przez nas wartości są małe i błąd przez to jest duży.

Teoria mówi, że sygnał oryginalny może zostać dokładnie zrekonstruowany ze swych próbek, jeśli próbki były pobierane w odstępach czasu T, zaś sygnał próbkowany nie zawiera częstotliwości większych niż połowa częstotliwości próbkowania. Czyli jeśli próbkujemy z częstotliwością 44100Hz, to sygnał nie powinien zawierać częstotliwości większych niż 22050Hz. W przeciwnym razie częstotliwości wyższe od limitu pojawią się w sygnale zrekonstruowanym zniekształcając go (aliasign). Twierdzenie Whittaker–Shannon'a mówi:

$x(t) = \sum_{n=-\infty}^{\infty} x[n] \cdot {\rm sinc}\left(\frac{t - nT}{T}\right)\,$

Inaczej sygnał oryginalny to splot impulsów z próbek z funkcją sinc o odpowiednim charakterze.

Teraz mając odtworzony sygnał ciągły, próbkujemy go ponownie z nową częstotliwością. Jeśli z zmniejszą to najpierw pozbywamy się z niego częstotliwości większych od połowy nowej częstotliwości próbkowania.

Całość robi poniższy kod. Jeśli mamy zrekonstruować sygnał z niższą częstotliwością to filtrowane są próbki wejściowe. Taki sygnał podlega rekonstrukcji zgodnie z wyżej podaną formułą tylko w interesujących nas czasach.

Kod funkcji:

function [ fout, tout ] = reconstruction( fin, fsin, fsout )

Tout = 1 / fsout;
Tin = 1 / fsin;

if (fsin == fsout)
    fout = fin;
    tout = Tout* linspace(0, size(fout, 2), size(fout, 2)); 
    return;
end

scale = fsout / fsin;

if (scale < 1)
    % for compatibility with downsampling and decimate
    fout = zeros(1, ceil(size(fin, 2)* scale));
else
    fout = zeros(1, fix(size(fin, 2) * scale));
end

tout = linspace(0, (1/fsin) * (size(fin, 2) - 1), size(fout, 2));            

% when downsampling cut frequencies below fsout/2 from fin 
% function.
if (fsin > fsout)
    k = sinc_gen(fsin, fsout/2);
    fin = conv(fin, k, 'same');
end

W = fix(160 / (fsout*2 / fsin)); % experimental
W = min(W, size(fin, 2));

tin = (0:1:size(fin, 2)-1) * Tin;

if (W > size(fin, 2))
    W = size(fin, 2);
end

for m=1:size(fout, 2)
    [~, min_index] = min(abs(tout(m) - tin));
    min_index = fix(min_index - W/2);
    if (min_index < 1)
        min_index = 1;
    end
    max_index = min_index + W - 1;
    if (max_index > size(fin, 2))
        max_index = size(fin, 2);
    end
    sc = sinc((tout(m) - (min_index-1:1:max_index-1)*Tin) / Tin);
    fout(m) = sum(fin(min_index:1:max_index) .* sc);
end

% for compatibility with upsampling
if (scale > 1)
    fout = fout(1:1:end-scale+1);
    tout = tout(1:1:end-scale+1);
end

Dla idealnej rekonstrukcji okno W powinno być równe długości funkcji, w praktyce może być znacznie krótsze. Jeśli mamy do czynienia z downsamplingiem to najpierw sygnał trzeba przefiltrować filtrem dolnoprzepustowym. Inaczej częstotliwości powyżej fsout, które były w fin odbiją się i pojawią się jako zakłócenia w sygnale. Z tak zakłóconym sygnałem już niż nie możemy zrobić - nie da się stwierdzić co jest oryginałem, a co zakłóceniem.

Rekonstrukcja taka jest dokładna, ale wyjątkowo czasochłonna, dla każdej próbki wyjściowej trzeba policzyć osobny impuls sinc.

Normalnie okno sinc powinno być równe wektorowi wejściowemu, tutaj jest krótsze, o długości powiedzmy eksperymentalnej dla której nie widać różnicy.

Jest też naniesionych parę poprawek na kompatybilność z innymi funkcjami dotyczącymi długości wektora wyjściowego.

Przykład:

%% input data

% resampling
fsin = 30000;
fsout = 7300;
Tin = 1/fsin;
scale = fsin / fsout;

func = @(ti) sin(2*pi*650*ti-pi/4) + sin(2*pi*1500*ti+pi/2);
tend = 1/500*2;

%build input func.
tin = (0:Tin:tend);
fin = func(tin);

% reconstruct
[fout, tout] = reconstruction(fin, fsin, fsout);

%resample
[fout1, tout1] = resampling(fin, fsin, fsout);

% calculate ideal reconstruction
forg = func(tout);

%% compare them

fig = figure;
hold on;
plot(tout, forg, 'r');
plot(0:Tin/20:tend, func(0:Tin/20:tend), 'g');
plot(tout, fout, 'b');
plot(tout1, fout1, 'black');
legend('ideal','in', 'out reconstr', 'out resampl');

%% compare fq

FORG = fftshift(abs(fft(forg)));
FIN = fftshift(abs(fft(func(0:Tin:tend))));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));

fqout = fsout/2*linspace(-1, 1, size(FORG, 2));
fqout1 = fsout/2*linspace(-1, 1, size(FOUT1, 2));
fqin = (1/Tin)/2*linspace(-1, 1, size(FIN, 2));

s = size(FIN, 2) / size(FOUT, 2);

hold on;
plot(fqout, s*FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, s*FOUT, 'b');
plot(fqout1, s*FOUT1, 'black');
legend('ideal','in', 'out reconstr', 'out resampl');

%% error

error_calc(fout, forg);
error_calc(fout1, forg);