close all
clear all
% Run the MATLAB filter design and analysis tool to design
% an FIR filter and export the coefficients to Num and Den variable

N = 2; % order of IIR filter
M = 3; % # of biquads

% IIR filter length is the filter order + 1
x = zeros(M, N+1); % input value (buffered)
y = zeros(M, N+1); % output values (buffered)
B = zeros(M, N+1); % numerator coefficients
A = zeros(M, N+1); % denominator coefficients
G = zeros(1,M+1); % scaling factors

fs = 48000;         % Sampling frequency
Ts = 1/fs;          % Sampling period
T = Ts:Ts:1;        % Time duration of signal

%% sum of 1 kHz cosine and its harmonics up to 24 kHz. 
f0 = 1000;
signal = cos(2*pi*f0*T);
for j = 2:24
    signal = signal + cos(2*pi*j*f0*T);
end

%% FIR Filter
y_filtered = zeros(1,length(T));   % Store output for one bufffer

for t = 1:length(T)
    x(1,1) = signal(t); % get input data sample as current input value
    
    for m = 1:M      
        %%% implementation of biquad(m)
    end

    y_filtered(t) = y(M,1);
    
end


%% Compute gain from input to output using peak-to-peak values
VppOutput = max(y_filtered(1+(N+1):end)) - min(y_filtered(1+(N+1):end));
VppInput= max(y_filtered(1:end-(N+1))) - min(y_filtered(1:end-(N+1)));
gain = VppOutput / VppInput;
gain_dB = 20*log10(gain); % gain in dB scale

figure
plot(T, y_filtered);
axis([-inf inf -1 1])

%% FFT
n = length(T);

figure
subplot(1,2,1)
f = fftshift(fft(signal));
plot(-fs/2:fs/n:(fs/2-fs/n), 10*log10(abs(f)))

subplot(1,2,2)
f = fftshift(fft(y_filtered));
plot(-fs/2:fs/n:(fs/2-fs/n), 10*log10(abs(f)))