%% Audio filters
%
% Table of Contents
% Basic lowpass and highpass filters
%     Signal time and spectrum plots
%     Signal low pass filtering
%     Signal High pass filtering
% Resonator (II order)
%     EXERCISE 1: implement the generalized resonator and plot the zero/pole patterns
%     POPCORN effect
% Equalizer and shelving filters
%     EXERCISE 2: implement a shelving LP filter with cut frew wc = 50 Hz (HP with wc = 4000 Hz )
%
%% Basic lowpass and highpass filters
%
%% Signal time and spectrum plots
% Plots ﻿in time show different graphical representations. The upper part shows about 8000 samples, the middle part the first 1000 samples and the lower part shows the first 100 samples out of a digital audio signal

% load file
clear;
filename = 'music/annoying.wav'
[y,fs] = audioread(filename);
info = audioinfo(filename)
plot(y)
xlabel('Sample Number')
ylabel('Amplitude')

% plot 1000 samples 
plot(y(5000:6000))
xlabel('Sample Number')
ylabel('Amplitude')

% plot 100 samples 
plot(y(10000:10100))
xlabel('Sample Number')
ylabel('Amplitude')

% play sound
sound(y,fs)

% spectrum 
Y = fft(y);
n = length(Y);           % number of samples
Py = abs(Y).^2/n;         % power of the DFT
Py = fftshift(Py);         % shift P values
f = (-n/2:n/2-1)*(fs/n); % 0-centered frequency range
plot(f,Py)
xlabel('Frequency')
ylabel('Power')

%
%% Signal low pass filtering
%
% The Matlab implementation for the filter is most easily accomplished using the filter function
% y = filter(B, A, x)
% The filter function takes three (3) arguments:
% feedforward coefficients B = [b0, a2, ..., bM]
% feedback coefficients A = [1, a1, a2 ,..., aN]
% the input signal x = [x0, x1, x2 ,..., xn]

% low pass Butterworth filter
fc = 2000;  % cut freq
[b, a] = butter(8, fc/(fs/2), 'low');  % normalized freq 
L = 2001;
[H,w] = freqz(b,a,L,fs,"whole");
plot(w-fs/2,fftshift(abs(H)))

% plot 
x = filter(b, a, y);
plot(x)
xlabel('Sample Number')
ylabel('Amplitude')

% play sound
sound(x,fs)

% spectrum
X = fft(x);
Px = abs(X).^2;        % power of the DFT
Px = fftshift(Px);         % shift y values
plot(f,Px)
xlabel('Frequency')
ylabel('Power')

%
%% Signal High pass filtering

% high pass Butterworth filter
fc = 10000;  % cut freq
[b, a] = butter(8, fc/(fs/2), 'high');  % normalized freq 
L = 2001;
[H,w] = freqz(b,a,L,fs,"whole");
plot(w-fs/2,fftshift(abs(H)))

% plot 
x = filter(b, a, y);
plot(x)
xlabel('Sample Number')
ylabel('Amplitude')

% play sound
sound(x,fs)

% spectrum
X = fft(x);
Px = abs(X).^2;        % power of the DFT
Px = fftshift(Px);         % shift y values
plot(f,Px)
xlabel('Frequency')
ylabel('Power')

%
%% Resonator (II order)
%
% Difference equation:
%
% Frequency response:
%
% Magnitude:
%
%
% Resonance frequency:
%

% main params
R = .99;
f0 = 10000; % resonance frequency

% angle of resonant phasor 
w0 = 2*pi*f0/fs;

% gain
G = (1-R)*sqrt(1-2*R*cos(2*w0)+R^2);
b = G;
a = [1 -2*R*cos(w0) R^2];

% freq response
freqz(b,a,L,fs)
[H,w] = freqz(b,a,L,fs);
plot(w,abs(H))
xlabel('Frequency')
ylabel('Power')
title('Linear scale')

% impulse response
impz(b,a,100,fs)

%
% Other filter example (see slides):

a = [1 -1.8831 0.9801];
b = 0.0062;
freqz(b,a,256,10000)
impz(b,a,160,10000)
zplane(b,a)

%
%% EXERCISE 1: implement the generalized resonator and plot the zero/pole patterns
%

%TODO...

%
%% POPCORN effect
% Single popcorn pop sound simulator: play a single pop with randomized properties
% Exponantially decaying white noise is filtered with a second-order resonator with center frequency about 400 Hz with some random variation and bandwidth between 50 and 150 Hz

 
% Filter a burst of noise x(n) using the resonator with zeros
% Bandwidth and center frequency are changedf ore very pop

fs = 20000;          % Sample rate (Hz)
fw = 10 + rand*50;  % Bandwidth of resonance (randomized 10...60 Hz)

fc = 100 + rand*250;  % frequency of resonance (randomized 50...300 Hz)

% Design resonator
B = 2*pi*fw/fs;             % Bandwidth in radians
R = 1 - (B/2);              % Radius of pole
a1 = -2*R*cos(2*pi*fc/fs);  % Coefficient a1
a2 = R^2;                   % Coefficient a2

% Design excitation signal
T = 10000;         % Duration in samples
exc = zeros(1,T);   % Initialize excitation sequence  
delta = [1 zeros(1,T-1)];  % delta signal d(n)
impulse = filter(1,[1 -0.999], delta);  % Exponential envelope
noiseseq = rand(1,T)-0.5;  % White noise sequence
exc(1:T) = impulse .* noiseseq;  % White noise with exponential envelope
plot(impulse)
xlabel('Time')
ylabel('Envelop')
stem(exc)

% spectrum of the envelop
Y = fft(exc);
n = length(Y);           % number of samples
Py = abs(Y).^2/n;         % power of the DFT
Py = fftshift(Py);         % shift P values
f = (-n/2:n/2-1)*(fs/n); % 0-centered frequency range
plot(f,Py)
xlabel('Frequency')
ylabel('Power')
title('spectrum of the impulse')

% Generate pop by filtering the noise sequence with the resonator
pop = filter([1 0 -1],[1 a1 a2],exc);
freqz([1 0 -1],[1 a1 a2],256,fs)
title('Resonator response')

plot(pop)
xlabel('Frequency')
ylabel('Power')
title('pop signal')

% spectrum 
Y = fft(pop);
n = length(Y);           % number of samples
Py = abs(Y).^2/n;         % power of the DFT
Py = fftshift(Py);         % shift P values
f = (-n/2:n/2-1)*(fs/n); % 0-centered frequency range
plot(f(n/2-1000:n/2+1000),Py(n/2-1000:n/2+1000))
xlabel('Frequency')
ylabel('Power')
title('spectrum of the pop')

% Play the pop
soundsc(pop,fs)

%
% Pop simulator: play a loop of pops with randomized properties

% Multiple popcorn pop sound simulator

 

for i = 1:20
    fw = 5 + rand*30;  
    fc = 50 + rand*100;  
    
    % Design resonator
    B = 2*pi*fw/fs;            
    R = 1 - (B/2);             
    a1 = -2*R*cos(2*pi*fc/fs); 
    a2 = R^2;                  
  
    impulse = filter(1,[1 -0.999],delta); 
    noiseseq = rand(1,T)-0.5;  
    exc(1:T) = impulse .* noiseseq; 
    
    % Generate pop by filtering the noise sequence with the resonator
    pop = filter([1 0 -1],[1 a1 a2],exc);
    
    soundsc(pop,fs)
    pause(rand(1))
end

%
%% Equalizer and shelving filters
%
% Computation of the coefficients for boost equalizer and  shelving filter with a given amplitude and cutoff frequency
% G    is the logrithmic gain (in dB)
% f0   is the center frequency
%  fs   is the sampling rate
% fw   bandwidth
% type define the filter type:  "EQ_boost", "EQ_cut", "SH_lp", "SH_hp"

% params
G0 = 1;     % reference
G = 9;      % dB
f0 = 200   % Hz
fs = 11025  % Hz
fw = 100    % Hz
type = "EQ_boost"   % "EQ_boost", "EQ_cut", "SH_lp", "SH_hp" 
% end params

G = 10^(G/20);      % gain in dB: 20*log(G/G0) = G dB
GB = 10^(-3/20)*G;  % gain in dB: 20*log(G/GB) = 3 dB
w0 = 2*pi*f0/fs;
Dw = 2*pi*fw/fs;
beta = sqrt((GB^2-G0^2)/(G^2-GB^2))*tan(Dw/2);
if type == "EQ_boost" || type == "EQ_cut"
    b0 = (G0+G*beta);
    b1 = -2*G0*cos(w0);
    b2 = (G0-G*beta);
    a1 = -2*cos(w0);
    a2 = (1-beta);
    b = [b0 b1 b2]/(1+beta);
elseif type == "SH_lp"
    b0 = (G0+G*beta);
    b1 = -2*G0*cos(w0);
    a1 = 2*G0*cos(w0);
    a2 = (1-beta);  
    b = [b0 b1]/(1+beta);
elseif type == "SH_hp"
    b0 = (G0+G*beta);
    b1 = -2*G0*cos(w0);
    a1 = 2*cos(w0);
    a2 = (1-beta);
    b = [b0 b1]/(1+beta);
end
a = [1 a1/(1+beta) a2/(1+beta)];
freqz(b,a,512,fs)
[H0,W0] = freqz(b,a,512,fs);
A0 = a;
B0 = b;

%
% Computation of the coefficients for cut equalizer and  shelving filter with a given amplitude and cutoff frequency
% G    is the logrithmic gain (in dB)
% f0   is the center frequency
%  fs   is the sampling rate
% fw   bandwidth
% type define the filter type:  "EQ_boost", "EQ_cut", "SH_lp", "SH_hp"

% params
G0 = 1;     % reference
G = -9;     % dB
f0 = 1000   % Hz
fs = 11025  % Hz
fw = 500    % Hz
type = "EQ_boost"   % "EQ_boost", "EQ_cut", "SH_lp", "SH_hp" 
% end params

G = 10^(G/20);      % gain in dB: 20*log(G/G0) = G dB
GB = 10^(3/20)*G;  % gain in dB: 20*log(G/GB) = 3 dB
w0 = 2*pi*f0/fs;
Dw = 2*pi*fw/fs;
beta = sqrt((GB^2-G0^2)/(G^2-GB^2))*tan(Dw/2);
if type == "EQ_boost" || type == "EQ_cut"
    b0 = (G0+G*beta);
    b1 = -2*G0*cos(w0);
    b2 = (G0-G*beta);
    a1 = -2*cos(w0);
    a2 = (1-beta);
    b = [b0 b1 b2]/(1+beta);
elseif type == "SH_lp"
    b0 = (G0+G*beta);
    b1 = -2*G0*cos(w0);
    a1 = 2*G0*cos(w0);
    a2 = (1-beta);  
    b = [b0 b1]/(1+beta);
elseif type == "SH_hp"
    b0 = (G0+G*beta);
    b1 = -2*G0*cos(w0);
    a1 = 2*cos(w0);
    a2 = (1-beta);
    b = [b0 b1]/(1+beta);
end
a = [1 a1/(1+beta) a2/(1+beta)];
freqz(b,a,512,fs);
[H1,W1] = freqz(b,a,512,fs);
A1 = a;
B1 = b;

%
% Plot the filter spectrum

plot(W0,20*log10(abs(H0)),'--g','LineWidth',2)
hold on
plot(W1,20*log10(abs(H1)),'--r','LineWidth',2)
plot(W1,20*log10(abs(H1))+20*log10(abs(H0)),'b','LineWidth',2) 
hold off

%
% Sound the music

filename = 'music/acoustic.wav'
[y,fs] = audioread(filename);
sound(y,fs)
% spectrum 
Y = fft(y);
n = length(Y);           % number of samples
Py = abs(Y).^2/n;         % power of the DFT
Py = fftshift(Py);         % shift P values
f = (-n/2:n/2-1)*(fs/n); % 0-centered frequency range
plot(f,Py)
xlabel('Frequency')
ylabel('Power')

%
% Apply the filter to a sound

x = filter(B0,A0,y);
x = filter(B1,A1,x);
sound(x,fs)
audiowrite("music/acoustic_EQ.wav",x,fs)
% spectrum 
Y = fft(x);
n = length(Y);           % number of samples
Py = abs(Y).^2/n;         % power of the DFT
Py = fftshift(Py);         % shift P values
f = (-n/2:n/2-1)*(fs/n); % 0-centered frequency range
plot(f,Py)
xlabel('Frequency')
ylabel('Power')

%
% MATLAB design tool for EQ filters

N = [2,2];
gain = [9,-9];
centerFreq = [0.05,0.25];
bwidth = [0.1,0.20];
[B,A] = designParamEQ(N,gain,centerFreq,bwidth);
SOS = [B',[ones(sum(N)/2,1),A']];
freqz(SOS,256,11025)

%
%% EXERCISE 2: implement a shelving LP filter with cut frew wc = 50 Hz (HP with wc = 4000 Hz )

% TODO...