Splot

Poniżej podane są funkcje w Matlabie realizujące splot, działające jak ich odpowiedniki w Matlabie:

function [ fout ] = conv_full( fin, k )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  fout = zeros(1,s1+s2-1);
  
  for m=1:size(fout, 2)
    fin1 = max(m-s2+1, 1);
    fin2 = min(m, s1);
    k1 = max(m-s1+1, 1);
    k2 = min(m, s2);
    fout(m) = sum(fin(fin1:fin2) .* k(k2:-1:k1));
  end

end

function [ fout ] = conv_same( fin, k )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  L = fix(size(k,2)/2)+1;
  H = L + size(fin, 2) - 1;
  
  fout = zeros(1, H-L+1);
  
  for m=L:H
    fin1 = max(m-s2+1, 1);
    fin2 = min(m, s1);
    k1 = max(m-s1+1, 1);
    k2 = min(m, s2);
    fout(m-L+1) = sum(fin(fin1:fin2) .* k(k2:-1:k1));
  end

end

Funkcja conv_same przydaje się w przetwarzaniu sygnałów, wtedy z reguły u traktujemy jako sygnał zaś v jako odwrotną transformatę Fouriera jakiegoś filtru, którego chcemy użyć. Splotowi u i v w czasie odpowiada mnożenie po stronie częstotliwości. Wynik conv_same będzie takiego samego rozmiaru jak u.

Funkcja v może być także odpowiedzią impulsową funkcji przejścia H(s) dowolnego układu. Znając H(s) liczymy jaka jest odpowiedź układu na impuls (1 po stronie s), czyli jest to transformata odwrotna 1*H(s). Mówimy wtedy o odpowiedzi impulsowej filtru. Po stronie czasu mamy układ o funkcji przejścia H(s), podajemy na jego wejście impuls (po stronie częstotliwości widmo impulsu jest ciągłe od minus nieskończoności do plus nieskończoności i o stałej amplitudzie). To co otrzymujemy na wyjściu to odpowiedź w czasie przefiltrowanych wszystkich możliwych częstotliwości. Transformatą Fouriera tej odpowiedzi wyznacza nam charakterystykę częstotliwością filtru.

I jeszcze przykład funkcji realizującej splot dla danego punktu:

function [ r ] = conv_point( fin, k, point )

  if (size(fin, 1) ~= 1)
    error('bad fin size');
  end
  
  if (size(k, 1) ~= 1)
    error('bad k size');
  end
  
  if (point < 1)
    error('bad point index');
  end
  
  if (point > size(fin, 2))
    error('bad point index');
  end
  
  s1 = size(fin, 2);
  s2 = size(k, 2);
  
  L = fix(size(k,2)/2)+1;
  H = L + size(fin, 2) - 1;
  
  m = point + L - 1;
  
  fin1 = max(m-s2+1, 1);
  fin2 = min(m, s1);
  k1 = max(m-s1+1, 1);
  k2 = min(m, s2);
  r = sum(fin(fin1:fin2) .* k(k2:-1:k1));

end

Porównanie wersji same i full:

fsin1 = 500; % czestotliwosc pierwsze sinusoidy
fsin1_phase = pi/2; %-pi/4.14; % przesuniecie fazowe pierwszej sinusoidy
fsin2 = 2000; % czestotliwosc drugiej sinusoidy
fsin2_phase = -pi/4; %pi/2.11; % przesuniecie fazowe drugiej sinusoidy
fs = 24000; % czestotliwosc probkowania

% funkcja wejsciowa.
func_in = @(ti) sin(2*pi*fsin1*ti+fsin1_phase) + ...
sin(2*pi*fsin2*ti+fsin2_phase);

fcut = 1000; % czestotliwosc odciecia
tend = 3/fsin1; % dlugosc funkcji
sinc_samples = 200;

% przygotuj wektor czasu.
ts = 1/fs; % odstep czasu pomiedzy probkowaniami
tx = (0:ts:tend);

if (mod(size(tx, 2), 2) == 1)
tx = tx(1:1:end-1);
end

% wejscie
fin = func_in(tx)/15;

% filtr
w1 = sinc_gen(fs, fcut, sinc_samples);

fout1 = conv(fin, w1, 'same');
fout2 = conv(fin, w1, 'full');

figure(1);
hold on;
plot(fin, 'r');
plot(w1, 'g');
plot(fout1, 'b');
plot(fout2, 'c');
legend('wejscie', 'filtr', 'conv full', 'conv same');

Przykład:

fs = 30000; % czestotliwosc probkowania
fcut = 1000; % czestotliwosc odciecia
fsin1 = 500; % czestotliwosc pierwsze sinusoidy
fsin1_phase = pi/2; %-pi/4.14; % przesuniecie fazowe pierwszej sinusoidy
fsin2 = 2000; % czestotliwosc drugiej sinusoidy
fsin2_phase = -pi/4; %pi/2.11; % przesuniecie fazowe drugiej sinusoidy
tend = 1/fsin1*3.35; % dlugosc funkcji
sinc_samples = 1000;

% funkcja wejsciowa.
func_in = @(ti) sin(2*pi*fsin1*ti+fsin1_phase) + ...
sin(2*pi*fsin2*ti+fsin2_phase);

% funkcja wyjsciowa - idealne odciecie.
func_out = @(ti) sin(2*pi*fsin1*ti+fsin1_phase);

% przygotuj wektor czasu.
ts = 1/fs; % odstep czasu pomiedzy probkowaniami
tx = (0:ts:tend);

% oblicz wartosci obu funkcji w czasie.
fin = func_in(tx);
fout = func_out(tx);

%[~, fin] = sinc_gen(fs, (fsin1+fsin2)*1.12, size(tx, 2));
%[~, fout] = sinc_gen(fs, fsin1, size(tx, 2));

% przygotuj sinca
[~, w1] = sinc_gen(fs, fcut, sinc_samples);

% oba okna w czestotliwosci.
fw1 = fft(w1);
fwx = linspace(-fs/2, fs/2, size(fw1, 2));

% funkcja wejsciowa w czestotliwosci.
ffin = fftshift(abs(fft(fin)));
ffinx = linspace(-fs/2, fs/2, size(ffin, 2));

% przyciecie poprzez splot
fout2 = conv(fin, w1, 'same');

figure(1);
hold on;
plot(fwx, fftshift(abs(fw1)), 'r');
legend('sinc');
title('filtr dolnoprzepustowy w przestrzeni czestotliwosci');

figure(2);
hold on;
plot(fin, 'c');
plot(fout);
plot(fout2, 'r');
title('porowanie w czasie');
legend('wejscie', 'idealne wyjscie', 'przyciecie przez splot');

figure(3);
hold on;
plot(ffinx, fftshift(abs(fft(fin))), 'c');
plot(ffinx, fftshift(abs(fft(fout))));
plot(ffinx, fftshift(abs(fft(fout2))), 'r');
title('porowanie w czestotliwosci');
legend('wejscie', 'idealne wyjscie', 'przyciecie przez splot');

Po odkomentowaniu linii zastępujących nasze funkcje sincami, wyniki przedstawiają się tak:

Test naszych funkcji:

% przetestujmy poprawnosc dzialania naszej funkcji conv_cut.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'same');
x12 = conv(a, b, 'same');
x13 = conv(a, c, 'same');
x14 = conv(a, d, 'same');

x21 = conv(b, a, 'same');
x22 = conv(b, b, 'same');
x23 = conv(b, c, 'same');
x24 = conv(b, d, 'same');

x31 = conv(c, a, 'same');
x32 = conv(c, b, 'same');
x33 = conv(c, c, 'same');
x34 = conv(c, d, 'same');

x41 = conv(d, a, 'same');
x42 = conv(d, b, 'same');
x43 = conv(d, c, 'same');
x44 = conv(d, d, 'same');

y11 = conv_same(a, a);
y12 = conv_same(a, b);
y13 = conv_same(a, c);
y14 = conv_same(a, d);

y21 = conv_same(b, a);
y22 = conv_same(b, b);
y23 = conv_same(b, c);
y24 = conv_same(b, d);

y31 = conv_same(c, a);
y32 = conv_same(c, b);
y33 = conv_same(c, c);
y34 = conv_same(c, d);

y41 = conv_same(d, a);
y42 = conv_same(d, b);
y43 = conv_same(d, c);
y44 = conv_same(d, d);

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

% przetestujmy poprawnosc dzialania naszej funkcji conv_full.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'full');
x12 = conv(a, b, 'full');
x13 = conv(a, c, 'full');
x14 = conv(a, d, 'full');

x21 = conv(b, a, 'full');
x22 = conv(b, b, 'full');
x23 = conv(b, c, 'full');
x24 = conv(b, d, 'full');

x31 = conv(c, a, 'full');
x32 = conv(c, b, 'full');
x33 = conv(c, c, 'full');
x34 = conv(c, d, 'full');

x41 = conv(d, a, 'full');
x42 = conv(d, b, 'full');
x43 = conv(d, c, 'full');
x44 = conv(d, d, 'full');

y11 = conv_full(a, a);
y12 = conv_full(a, b);
y13 = conv_full(a, c);
y14 = conv_full(a, d);

y21 = conv_full(b, a);
y22 = conv_full(b, b);
y23 = conv_full(b, c);
y24 = conv_full(b, d);

y31 = conv_full(c, a);
y32 = conv_full(c, b);
y33 = conv_full(c, c);
y34 = conv_full(c, d);

y41 = conv_full(d, a);
y42 = conv_full(d, b);
y43 = conv_full(d, c);
y44 = conv_full(d, d);

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

%% przetestujmy poprawnosc dzialania naszej funkcji conv_point.

a = [1,2,3,4];
b = [1,2,3,4,5];
c = [1,2,3,4,5,6];
d = [1,2,3,4,5,6,7];

x11 = conv(a, a, 'same');
x12 = conv(a, b, 'same');
x13 = conv(a, c, 'same');
x14 = conv(a, d, 'same');

x21 = conv(b, a, 'same');
x22 = conv(b, b, 'same');
x23 = conv(b, c, 'same');
x24 = conv(b, d, 'same');

x31 = conv(c, a, 'same');
x32 = conv(c, b, 'same');
x33 = conv(c, c, 'same');
x34 = conv(c, d, 'same');

x41 = conv(d, a, 'same');
x42 = conv(d, b, 'same');
x43 = conv(d, c, 'same');
x44 = conv(d, d, 'same');

y11 = zeros(1, size(a, 2));
y12 = zeros(1, size(a, 2));
y13 = zeros(1, size(a, 2));
y14 = zeros(1, size(a, 2));
for i=1:1:size(a, 2)
y11(i) = conv_point(a, a, i);
y12(i) = conv_point(a, b, i);
y13(i) = conv_point(a, c, i);
y14(i) = conv_point(a, d, i);
end

y21 = zeros(1, size(b, 2));
y22 = zeros(1, size(b, 2));
y23 = zeros(1, size(b, 2));
y24 = zeros(1, size(b, 2));
for i=1:1:size(b, 2)
y21(i) = conv_point(b, a, i);
y22(i) = conv_point(b, b, i);
y23(i) = conv_point(b, c, i);
y24(i) = conv_point(b, d, i);
end

y31 = zeros(1, size(c, 2));
y32 = zeros(1, size(c, 2));
y33 = zeros(1, size(c, 2));
y34 = zeros(1, size(c, 2));
for i=1:1:size(c, 2)
y31(i) = conv_point(c, a, i);
y32(i) = conv_point(c, b, i);
y33(i) = conv_point(c, c, i);
y34(i) = conv_point(c, d, i);
end

y41 = zeros(1, size(d, 2));
y42 = zeros(1, size(d, 2));
y43 = zeros(1, size(d, 2));
y44 = zeros(1, size(d, 2));
for i=1:1:size(d, 2)
y41(i) = conv_point(d, a, i);
y42(i) = conv_point(d, b, i);
y43(i) = conv_point(d, c, i);
y44(i) = conv_point(d, d, i);
end

test = @(a,b) assert((size(a,2) == size(b,2)) && (sum(abs(a - b)) == 0));

test(x11, y11);
test(x12, y12);
test(x13, y13);
test(x14, y14);

test(x21, y21);
test(x22, y22);
test(x23, y23);
test(x24, y24);

test(x31, y31);
test(x32, y32);
test(x33, y33);
test(x34, y34);

test(x41, y41);
test(x42, y42);
test(x43, y43);
test(x44, y44);

Funkcja pomocnicza:

function [ k, t ] = sinc_gen( sampling_rate, cutfq, samples )

T = 1 / sampling_rate;

if (nargin == 2)
samples = 160 / (cutfq*2 / sampling_rate); % experimental
end

tout = T * (1:samples);
tout = tout - T*(samples+1)/2;
B = 2*cutfq;
k = sinc(B*tout);
k = k / trapz(k);

if (nargout == 2)
t = tout;
end

end