function y = ka_FIR(x, fs, M)
% KA_FIR keskiarvoistaa M:llä pisteellä
%
% y = ka_FIR(x, fs, M)
%     y    suodatettu sekvenssi
%     x    suodatettava sekvenssi
%     fs   näytteenottotaajuus
%     M    liikkuva M:n pisteen keskiarvoistava suodin
%
% Esimerkki:
%     [x, fs] = wavread('aani.wav');
%     ka_FIR(x, fs, 3);  % kolmen pisteen suodin

% Jukka Parviainen, 2004, 21.2.2005

figure(1); clf;
t = [1 : length(x)]/fs;
y = zeros(1, length(x));
%
% Voidaan tehdä silmukassa:
%
%for k = M:length(x)
%  y(k) = (1/M) * sum(x(k-M+1 : k));
%end;
%
% mutta *sama* Matlabissa tehokkaammin toteutettu:
%
y = filter((1/M)*ones(1,M), 1, x);
%
%
figure(1), clf;
plot(t, x, 'b', t, y, 'r');
grid on;
xlabel('Aika (s)');
legend({'Alkuper.','Suodatettu'});
title('Aikataso');

figure(2); clf;
xF= fft(x(45000:44999+2^10));  % vain pieni pala
maxx = max(abs(xF(:)));
xF= xF/maxx;
yF= fft(y(45000:44999+2^10))/maxx;  % vain pieni pala
f = fs*[0 : 1023]/1024;
plot(f, abs(xF), 'b');
hold on;
plot(f, abs(yF), 'r');

w = [0 : pi/2048 : pi];  % omega/rad-asteikko
W = (fs/2)*w/pi;         % taajuus/Hz-asteikko
H = (1/M)*sum(exp(-j*[0:M-1]'*w));
plot(W, abs(H), 'k');         % amplitudivaste 
grid on;
xlabel('Taajuus (Hz)');
axis([0 fs/2 -0.1 1.2]);
title('Taajuustaso');
legend({'Alkuper.','Suodatettu','Suodin'});
