TP3_8: Noise removal done.

TP3/DiogoEliseu_TP3_7.m

@@ -10,6 +10,7 @@ s = y(:,1);
 % F: vector of cyclical frequencies
 % T: vector of time instants
 % P: power spectral density (PSD)
+                    % signal, window, noverlap, nfft, fs
 [S, F, T, P] = spectrogram(s, 128, 96, 128, Fs);
 
 figure(1)

TP3/DiogoEliseu_TP3_8.m

@@ -1,34 +0,0 @@
-%% clean environment
-clear; close all; clc
-
-%% read input
-[signal, fs] = audioread("Canto1.mp3");
-
-%% Compute the FFT
-signal_length = length(signal);
-fhat = fft(signal, signal_length);
-PSD = fhat.*conj(fhat)/signal_length; % Power spectrum density
-%[S, freq, T, PSDlel] = spectrogram(signal, 128, 96, 128, Fs);
-% signal, window, noverlap, nfft, fs
-freq_vector = linspace(1, fs, signal_length);
-L = 1:floor(signal_length/2);
-
-cutoff_indices = abs(PSD)>0.002;
-PSD_clean = PSD.*cutoff_indices;
-fhat_filtered = cutoff_indices.*fhat;
-signal_filtered = ifft(fhat_filtered);
-
-figure(1)
-plot(PSD)
-
-figure(2)
-plot(PSD_clean)
-
-figure(3)
-plot(signal)
-
-figure(4)
-plot(signal_filtered)
-
-%% output
-audiowrite('restored.flac', signal_filtered, fs);

TP3_8/DiogoEliseu_TP3_8_butterworth.m

@@ -0,0 +1,58 @@
+%% Inicialização do ambiente
+clear; close all; clc
+
+%% >Exercício 8<
+%% Input
+[y, Fs] = audioread("Canto1.mp3"); % Signal and Sampling Frequency
+signal = y(:,1);                   % Discard second channel
+Fn = Fs/2;                         % Nyquist Frequency (Hz)
+L = length(signal);                % Signal Length
+
+%% Find the noise
+dft_signal = fft(signal)./L;          % Fourier Transform
+Fv = linspace(0, 1, fix(L/2)+1)*Fn;   % Frequency Vector
+
+%% Filter
+%signal_filtered = filter(DiogoEliseu_TP3_8_butterworth_2nd_order_filter,signal);
+signal_filtered = filter(DiogoEliseu_TP3_8_butterworth_filter,signal);
+
+%% Filtered signal
+dft_signal_filtered = fft(signal_filtered)/L; % Fourier Transform
+
+%% Output
+signal_power = abs(dft_signal(1:(L/2+1)));
+figure(1)
+subplot(2,3,1)
+hold on
+plot(Fv, signal_power)
+grid
+marks = [2850, 3020, 7000, 7150];
+plot(2850, signal_power(2850), 'r*')
+plot(3020, signal_power(3020), 'r*')
+plot(7000, signal_power(7000), 'r*')
+plot(7150, signal_power(7150), 'r*')
+xlabel("Frequência (Hz)"); ylabel("|H_{1}(f)|"); title("DFT do Sinal")
+hold off
+subplot(2,3,4)
+plot(Fv, abs(dft_signal_filtered(1:(L/2+1))))
+grid
+xlabel("Frequência (Hz)"); ylabel("|H_{1}(f)|"); title("DFT do Sinal Filtrado")
+subplot(2,3,2)
+plot(Fv, signal(1:(L/2+1)))
+grid
+xlabel("Frequência (Hz)"); ylabel("Amplitude"); title("Sinal")
+subplot(2,3,5)
+plot(Fv, signal_filtered(1:(L/2+1)))
+grid
+xlabel("Frequência (Hz)"); ylabel("Amplitude"); title("Sinal Filtrado")
+subplot(2,3,3)
+plot(Fv, mag2db(signal_power))
+grid
+xlabel("Frequência (Hz)"); ylabel("Ganho (dB)"); title("Ganho do Sinal")
+subplot(2,3,6)
+plot(Fv, mag2db(abs(dft_signal_filtered(1:(L/2+1)))))
+grid
+xlabel("Frequência (Hz)"); ylabel("Ganho (dB)"); title("Ganho do Sinal Filtrado")
+
+%sound(signal_filtered, Fs);
+%audiowrite('canto1_butterworth.flac', signal_filtered, Fs);

TP3_8/DiogoEliseu_TP3_8_butterworth_2nd_order_filter.m

@@ -0,0 +1,21 @@
+function Hd = DiogoEliseu_TP3_8_final2
+%DIOGOELISEU_TP3_8_FINAL2 Returns a discrete-time filter object.
+
+% MATLAB Code
+% Generated by MATLAB(R) 9.8 and Signal Processing Toolbox 8.4.
+% Generated on: 31-Dec-2020 16:39:59
+
+% Butterworth Bandpass filter designed using FDESIGN.BANDPASS.
+
+% All frequency values are in Hz.
+Fs = 44100;  % Sampling Frequency
+
+N   = 2;     % Order
+Fc1 = 2850;  % First Cutoff Frequency
+Fc2 = 7150;  % Second Cutoff Frequency
+
+% Construct an FDESIGN object and call its BUTTER method.
+h  = fdesign.bandpass('N,F3dB1,F3dB2', N, Fc1, Fc2, Fs);
+Hd = design(h, 'butter');
+
+% [EOF]

TP3_8/DiogoEliseu_TP3_8_butterworth_filter.m

@@ -0,0 +1,27 @@
+function Hd = DiogoEliseu_TP3_8_final_filter
+%DIOGOELISEU_TP3_8_FINAL_FILTER Returns a discrete-time filter object.
+
+% MATLAB Code
+% Generated by MATLAB(R) 9.8 and Signal Processing Toolbox 8.4.
+% Generated on: 31-Dec-2020 19:24:38
+
+% Butterworth Bandpass filter designed using FDESIGN.BANDPASS.
+
+% All frequency values are in Hz.
+Fs = 44100;  % Sampling Frequency
+
+Fstop1 = 2850;        % First Stopband Frequency
+Fpass1 = 3020;        % First Passband Frequency
+Fpass2 = 7000;        % Second Passband Frequency
+Fstop2 = 7150;        % Second Stopband Frequency
+Astop1 = 96;          % First Stopband Attenuation (dB)
+Apass  = 1;           % Passband Ripple (dB)
+Astop2 = 104;         % Second Stopband Attenuation (dB)
+match  = 'stopband';  % Band to match exactly
+
+% Construct an FDESIGN object and call its BUTTER method.
+h  = fdesign.bandpass(Fstop1, Fpass1, Fpass2, Fstop2, Astop1, Apass, ...
+                      Astop2, Fs);
+Hd = design(h, 'butter', 'MatchExactly', match);
+
+% [EOF]

TP3_8/DiogoEliseu_TP3_8_high_pass.m

@@ -0,0 +1,55 @@
+%% Inicialização do ambiente
+clear; close all; clc
+
+%% >Exercício 8<
+
+%% Input
+[y, Fs] = audioread("Canto1.mp3"); % Signal and Sampling Frequency
+signal = y(:,1);                   % Discard second channel
+Fn = Fs/2;                         % Nyquist Frequency (Hz)
+L = length(signal);                % Signal Length
+
+%% Compute the DFT
+dft = fft(signal)./L;
+magnitude = 10.*log10(abs(dft).^2); % power is squared, thus 10
+Fv = linspace(0, 1, L)*Fn;          % Frequency Vector
+
+%% High-pass Filter
+not_to_cutoff_indices = magnitude > -96;
+dft_filtered = not_to_cutoff_indices.*dft;
+magnitude_filtered = not_to_cutoff_indices.*magnitude;
+signal_filtered = ifft(dft_filtered.*L);
+
+%% Output
+Fv_half = linspace(0, 1, fix(L/2)+1)*Fn;
+y = abs(dft(1:(L/2+1)));
+y_filtered = abs(dft_filtered(1:(L/2+1)));
+figure(1)
+subplot(2,3,1)
+plot(Fv_half, y)
+grid
+xlabel("Frequência (Hz)"); ylabel("|H_{1}(f)|"); title("DFT do Sinal")
+subplot(2,3,4)
+plot(Fv_half, y_filtered)
+grid
+xlabel("Frequência (Hz)"); ylabel("|H_{1}(f)|"); title("DFT do Sinal Filtrado")
+subplot(2,3,2)
+plot(Fv_half, signal(1:(L/2+1)))
+grid
+xlabel("Frequência (Hz)"); ylabel("Amplitude"); title("Sinal")
+subplot(2,3,5)
+plot(Fv_half, y_filtered)
+grid
+xlabel("Frequência (Hz)"); ylabel("Amplitude"); title("Sinal Filtrado")
+subplot(2,3,3)
+%plot(Fv, magnitude);
+plot(Fv_half, mag2db(y))
+grid
+xlabel("Frequência (Hz)"); ylabel("Ganho (dB)"); title("Ganho do Sinal")
+subplot(2,3,6)
+plot(Fv_half, mag2db(y_filtered))
+grid
+xlabel("Frequência (Hz)"); ylabel("Ganho (dB)"); title("Ganho do Sinal Filtrado")
+
+%sound(signal_filtered, Fs);
+%audiowrite('canto1_high_pass.flac', signal_filtered, Fs);