clear all
close all 
clc

% Modify the desired frequency for the whole code here
desired_freq = 4000;
% Read and Set up coeficient
[y, fs] = audioread("sound/modulator22.wav");
decrease_coef = round(fs/desired_freq);

% Modify Signal with downsample
% Only keeps one sample out of decrease_coef
down_y = downsample(y,decrease_coef);
down_fs = fs/decrease_coef;
audiowrite("sound/down_output.wav", down_y, down_fs);
% Plot Down Modified 
plot(0:1/down_fs:(length(down_y)-1)/down_fs,down_y);
xlabel("Time (s)");
ylabel("Amplitude");
title("Downsample Sound Amplitude Over Time");

% Modify Signal with decimate
% Only keeps one sample out of decrease_coef
deci_y = decimate(y,decrease_coef);
deci_fs = fs/decrease_coef;
audiowrite("sound/deci_output.wav", deci_y, deci_fs);
% Plot Down Modified 
figure;
plot(0:1/deci_fs:(length(deci_y)-1)/deci_fs,deci_y);
xlabel("Time (s)");
ylabel("Amplitude");
title("Decimate Sound Amplitude Over Time");


% Filter Signal with FIR filter
order = 30;
cutoff_freq = 1000; % in Hz
fir_b = fir1(order, cutoff_freq/(fs/2));
% Show Filter
figure
freqz(fir_b)
% Apply filter
fir_y = filter(fir_b, 1, y);
audiowrite("sound/fir_output.wav", fir_y, fs);
% Plot Down Modified 
figure;
plot(0:1/fs:(length(fir_y)-1)/fs,fir_y);
xlabel("Time (s)");
ylabel("Amplitude");
title("FIR Filter Sound Amplitude Over Time");

% Filter Signal with Butter filter
order = 8; 
cutoff_freq = 1000; % in Hz
[b, a] = butter(order, cutoff_freq/(fs/2), 'low');
% Show Filter
figure
freqz(b)
% Apply filter
butt_y = filter(b, a, y);
audiowrite("sound/butt_output.wav", butt_y, fs);
% Plot Down Modified 
figure;
plot(0:1/fs:(length(butt_y)-1)/fs,butt_y);
xlabel("Time (s)");
ylabel("Amplitude");
title("Butter Filter Sound Amplitude Over Time");

% Test FIR Filter Stability
% Find the transfer function
H = tf(fir_b, 1);
% Get the poles of the transfer function
poles = pole(H);
% Check if all poles are inside the unit circle
if all(abs(poles) < 1)
    disp('The FIR filter is stable');
else
    disp('The FIR filter is unstable');
end
% Test Butter Filter Stability
% Find the transfer function
H = tf(b, a);
% Get the poles of the transfer function
poles = pole(H);
% Check if all poles are inside the unit circle
if all(abs(poles) < 1)
    disp('The Butter filter is stable');
else
    disp('The Butter filter is unstable');
end

% Modify Signal with downsample
% Only keeps one sample out of decrease_coef
down_fir_y = downsample(fir_y,decrease_coef);
down_fir_fs = fs/decrease_coef;
audiowrite("sound/down_fir_output.wav", down_fir_y, down_fir_fs);
% Plot Down Modified 
figure
plot(0:1/down_fir_fs:(length(down_fir_y)-1)/down_fir_fs,down_fir_y);
xlabel("Time (s)");
ylabel("Amplitude");
title("Downsample FIR Sound Amplitude Over Time");

% Modify Signal with downsample
% Only keeps one sample out of decrease_coef
down_butt_y = downsample(butt_y,decrease_coef);
down_butt_fs = fs/decrease_coef;
audiowrite("sound/down_butt_output.wav", down_butt_y, down_butt_fs);
% Plot Down Modified 
figure
plot(0:1/down_butt_fs:(length(down_butt_y)-1)/down_butt_fs,down_butt_y);
xlabel("Time (s)");
ylabel("Amplitude");
title("Downsample Butter Sound Amplitude Over Time");