%   This code simulates a running average of 16 data points.   Download and run.

Fs=10e3;        %   10 kHz sampling frequency

T=1/Fs;
N=4096;          %   N samples  in a Unit Test Period 
M=16;					%   Average 16 samples
freq = Fs/N:Fs/N:Fs;       %   Frequencies are multiples of the fundamental frequency Ff = Fs/N
w = 2*pi*freq;

z=exp(j*w*T);           %   Z-Transform variable 

hz = (1/M)*(1-z.^-M)./(1-z.^-1);   %   Transfer function where the output is 
                                   %   a running average of 16 input samples
                                   %    (Theory)

%   Set up a multitone signal to test the averager.
testbins =[2 4 10 20 22 24 25 26 27 28 30 40 60 120];
testbins=10*testbins;       %  Adjust test frequencies.  Test freq = testbins(k)*Fs/N
%   Test the filter with length(testbins)=14 tones 
x=zeros(1,N);         %   Initilize input tones, x at zero

for k=1:length(testbins)
    for n=1:N
    x(n) = x(n) + cos(pi/2 + 2*pi*testbins(k)*(n-1)/N); %   Tones have 90 degree phase so x(1)=0
    end
end

y=zeros(1,N);

for n=1:N
   if n<=M
      for k=0:n-1
         y(n)=y(n)+x(n-k)/n;
      end
   else     
   for k=1:M
      y(n) = y(n)+(1/M)*x(n-k);   %   Averager filter
   end
end
end

L=N/2;

yf=2*fft(y)/length(y);   %   Take the Fourier transform of y(n)
yf(1)=yf(1);             %   Correct the DC term. The term at freq = Fs/2 should also be corrected.
ym = abs(yf);            %   Get amplitude
ydB=20*log10(ym);        %   Convert to dB

semilogx(freq(1:L),abs(hz(1:L)),'-',freq(testbins+1),ym(testbins+1),'o') % Plot 
ylabel('Filter gain (Vout/Vin)'), xlabel('frequency (Hz)')
legend('Theory','Measured')  
%  The magnitude of hz is the theory. ym, the magnitude of yf. is the measured.

