pkg load signal

% Load the audio signal
[signal, Fs] = audioread("Sound/modulator22.wav");
t = (0:length(signal)-1) / Fs; % Time vector

% Design low-pass filter
lpf_order = 30;
lpf_cutoff_freq = 1000; % in Hz
lpf_filter = fir1(lpf_order, lpf_cutoff_freq / (Fs/2));

% Apply low-pass filter to the original signal
filtered_signal_lpf = filter(lpf_filter, 1, signal);

% Downsample the filtered signal using downsample()
downsampled_signal1 = downsample(filtered_signal_lpf, round(Fs/4000));

% Downsample the filtered signal using decimate()
downsampled_signal2 = decimate(filtered_signal_lpf, round(Fs/4000));

% Write downsampled signals to audio files
audiowrite('downsampled_signal1.wav', downsampled_signal1, 4000);
audiowrite('downsampled_signal2.wav', downsampled_signal2, 4000);

% Parameters for measurements
num_measurements = 100; % Number of measurements
durations_dft = zeros(1, num_measurements);
durations_fft = zeros(1, num_measurements);

for i = 1:num_measurements
    % Measure time taken for DFT
    tic;
    [power_dft, duration_dft] = frequencySpectrum_old(signal, Fs, false, false);
    durations_dft(i) = duration_dft;

    % Measure time taken for FFT with padding
    tic;
    [power_fft, duration_fft] = frequencySpectrum_old(signal, Fs, true, false);
    durations_fft(i) = duration_fft;
end

% Calculate average and standard deviation
avg_duration_dft = mean(durations_dft);
std_dev_dft = std(durations_dft);

avg_duration_fft = mean(durations_fft);
std_dev_fft = std(durations_fft);

fprintf('Average duration for DFT: %f seconds\n', avg_duration_dft);
fprintf('Standard deviation for DFT: %f seconds\n', std_dev_dft);
fprintf('\n');
fprintf('Average duration for FFT (with padding): %f seconds\n', avg_duration_fft);
fprintf('Standard deviation for FFT (with padding): %f seconds\n', std_dev_fft);

% Apply FIR filter
% Design FIR filter
fir_order = 30;
cutoff_freq_fir = 1000; % in Hz
fir_filter = fir1(fir_order, cutoff_freq_fir / (Fs/2));

% Apply FIR filter to the original signal
filtered_signal_fir = filter(fir_filter, 1, signal);

% Check FIR filter stability (FIR filters are always stable)

% Apply IIR filter
% Design Butterworth IIR filter
iir_order = 8;
cutoff_freq_iir = 1000; % in Hz
[b, a] = butter(iir_order, cutoff_freq_iir / (Fs/2), 'low');

% Check IIR filter stability
poles = roots(a); % Compute poles
iir_stable = all(abs(roots(a)) < 1);  % Check if all poles are inside the unit circle
fprintf('IIR filter stability: %d\n', iir_stable);
fprintf('Poles of the IIR filter: ');
disp(poles');

% Apply IIR filter to the original signal
filtered_signal_iir = filter(b, a, signal);

% Plot all signals and filtered signals on the same figure
figure;

% Original Signal
subplot(4, 1, 1);
plot(t, signal);
xlabel('Time (s)');
ylabel('Amplitude');
title('Original Signal');

% Low-pass Filtered Signal
subplot(4, 1, 2);
plot(t, filtered_signal_lpf);
xlabel('Time (s)');
ylabel('Amplitude');
title('Low-pass Filtered Signal');

% Downsampled Signal using downsample()
subplot(4, 1, 3);
t_down1 = (0:length(downsampled_signal1)-1) / (Fs/4000);
plot(t_down1, downsampled_signal1);
xlabel('Time (s)');
ylabel('Amplitude');
title('Downsampled Signal using downsample()');

% Downsampled Signal using decimate()
subplot(4, 1, 4);
t_down2 = (0:length(downsampled_signal2)-1) / (Fs/4000);
plot(t_down2, downsampled_signal2);
xlabel('Time (s)');
ylabel('Amplitude');
title('Downsampled Signal using decimate()');


% Plot filtered signals
figure;

% FIR Filtered Signal
subplot(2, 1, 1);
plot(t, filtered_signal_fir);
xlabel('Time (s)');
ylabel('Amplitude');
title('FIR Filtered Signal');

% IIR Filtered Signal
subplot(2, 1, 2);
plot(t, filtered_signal_iir);
xlabel('Time (s)');
ylabel('Amplitude');
title('IIR Filtered Signal');

% Display frequency response of FIR and IIR filters
figure;

% Frequency response of FIR filter
[h_fir, w_fir] = freqz(fir_filter, 1, 1024, Fs);
subplot(2, 1, 1);
plot(w_fir, 20*log10(abs(h_fir)));
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Frequency Response of FIR Filter');
grid on;

subplot(2, 1, 2);
plot(w_fir, angle(h_fir));
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');
title('Phase Response of FIR Filter');
grid on;

% Frequency response of IIR filter
[h_iir, w_iir] = freqz(b, a, 1024, Fs);
figure;
subplot(2, 1, 1);
plot(w_iir, 20*log10(abs(h_iir)));
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Frequency Response of IIR Filter');
grid on;

subplot(2, 1, 2);
plot(w_iir, angle(h_iir));
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');
title('Phase Response of IIR Filter');
grid on;

% Downsample filtered signals to 4000 Hz
downsampled_filtered_signal_fir = downsample(filtered_signal_fir, round(Fs/4000));
downsampled_filtered_signal_iir = downsample(filtered_signal_iir, round(Fs/4000));

% Compare downsampled filtered signals
figure;

% FIR Filtered Signal (downsampled)
subplot(2, 1, 1);
plot(t_down1, downsampled_filtered_signal_fir);
xlabel('Time (s)');
ylabel('Amplitude');
title('FIR Filtered Signal (Downsampled)');

% IIR Filtered Signal (downsampled)
subplot(2, 1, 2);
plot(t_down1, downsampled_filtered_signal_iir);
xlabel('Time (s)');
ylabel('Amplitude');
title('IIR Filtered Signal (Downsampled)');

% Write downsampled signals to audio files
audiowrite('downsampled_signal1_fir.wav', downsampled_filtered_signal_fir, 4000);
audiowrite('downsampled_signal2_iir.wav', downsampled_filtered_signal_iir, 4000);