Upsampling

Upsampling czyli zwiększenie częstotliwości próbkowania sygnału o zadany całkowity współczynik przy zachowaniu jego widma.

Funkcja:

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

scale = fsout / fsin;

if (fix(scale) ~= scale)
    error('scale should be integer');
end

if (scale <= 0)
    error('scale should be greater than zero');
end

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

[~, k] = sinc_gen(fsout, fsin/2);
if (mod(size(k, 2), 2) == 0)
    [~, k] = sinc_gen(fsout, fsin/2, fix(size(k, 2) / 2) * 2 + 1);
end

fout = conv(fout, k, 'same');

fout = fout(1:1:end-scale+1);
tout = tout(1:1:end-scale+1);

end

Funkcja upsample uzupełnia fin zerami. Zera w ilości scale-1 są wstawiane między próbki i na koniec. My te końcowe zera usuwamy. Dzięki czemu ostatnia próbka sygnału wyjściowego odpowiada ostatniej próbce sygnału wejściowego.

Okno sinc musi mieć nieparzystą długość by uniknącć przesunięcia wyniku.

Sam splot to tzw. idealna rekonstrukcja, tutaj o tyle prostsza i szybsza, że współczynnik skalowania jest całkowity i sinc rekonstruujący może być wyliczony raz i używamy dla każdego punktu z osobna.

W stosunku do wbudowanej w Matlab funkcji interp błąd przez nas popełniany jest na tyle niewielki, że możemy powiedzieć, że wynika on z innej metody interpolacji.

Przykład:

%% input data

fsin = 6000;
fsout = 30000;
scale = fsout / fsin;

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

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

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

% calculate ideal reconstruction
forg = func(tout);

fout1 = interp(fin, scale);
fout1 = fout1(1:1:end-scale+1);

%% compare them

hold on;
plot(tout, forg, 'r');
plot(tin, fin, 'g');
plot(tout, fout, 'b');
plot(tout, fout1, 'm');
legend('ideal','in', 'out', 'interp');

%% compare fq

FIN = fftshift(abs(fft(fin)));
FOUT = fftshift(abs(fft(fout)));
FOUT1 = fftshift(abs(fft(fout1)));
s = size(FIN, 2) / size(FOUT, 2);
FORG = s*fftshift(abs(fft(forg)));
FOUT = s*FOUT;
FOUT1 = s*FOUT1;

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

hold on;
plot(fqout, FORG, 'r');
plot(fqin, FIN, 'g');
plot(fqout, FOUT, 'b');
plot(fqout, FOUT1, 'm');
legend('ideal','in', 'out', 'interp');

%% error

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