That's it for noise removal

TP3_8/DiogoEliseu_TP3_8_high_pass.m

@@ -10,44 +10,45 @@ 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
+dft_signal = fft(signal)./L;
+
+%% Compute the Magninude in dB from Power
+magnitude = 10.*log10(abs(dft_signal).^2); % power is squared, thus 10
 
 %% High-pass Filter
 not_to_cutoff_indices = magnitude > -96;
-dft_filtered = not_to_cutoff_indices.*dft;
+dft_filtered = not_to_cutoff_indices.*dft_signal;
 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)));
+Fv = linspace(0, 1, fix(L/2)+1)*Fn;
+y = abs(dft_signal(1:(L/2+1)));
 y_filtered = abs(dft_filtered(1:(L/2+1)));
 figure(1)
 subplot(2,3,1)
-plot(Fv_half, y)
+plot(Fv, y)
 grid
 xlabel("Frequência (Hz)"); ylabel("|H_{1}(f)|"); title("DFT do Sinal")
 subplot(2,3,4)
-plot(Fv_half, y_filtered)
+plot(Fv, 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)))
+plot(Fv, signal(1:(L/2+1)))
 grid
 xlabel("Frequência (Hz)"); ylabel("Amplitude"); title("Sinal")
 subplot(2,3,5)
-plot(Fv_half, signal_filtered(1:(L/2+1)))
+plot(Fv, signal_filtered(1:(L/2+1)))
 grid
 xlabel("Frequência (Hz)"); ylabel("Amplitude"); title("Sinal Filtrado")
 subplot(2,3,3)
 %plot(Fv, magnitude);
-plot(Fv_half, mag2db(y))
+plot(Fv, mag2db(y))
 grid
 xlabel("Frequência (Hz)"); ylabel("Ganho (dB)"); title("Ganho do Sinal")
 subplot(2,3,6)
-plot(Fv_half, mag2db(y_filtered))
+plot(Fv, mag2db(y_filtered))
 grid
 xlabel("Frequência (Hz)"); ylabel("Ganho (dB)"); title("Ganho do Sinal Filtrado")
 

TP3_8/trash/DiogoEliseu_TP3_8_butter.m

@@ -0,0 +1,50 @@
+%% Inicialização do ambiente
+clear; close all; clc
+
+%% Exercício 8
+% Input
+[y, Fs] = audioread("Canto1.mp3"); % Signal and Sampling Frequency
+signal = y(:,1);
+Fn = Fs/2; % Nyquist Frequency (Hz)
+%Ts = 1/Fs;
+L = length(signal); % Signal Length
+
+% Find the noise
+FTsignal = fft(signal)/L;           % Fourier Transform
+Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector
+Iv = 1:numel(Fv);                   % Index Vector
+figure(1)
+% x is Hz, y is amplitude
+plot(Fv, abs(FTsignal(Iv))*2)
+grid
+
+% Butterworth Low-Pass Order 7
+%[num, den] = butter(7, fc, 'low'); %order, cutoff frequency, type
+%[SOS,G] = tf2sos(num, den); % Convert To Second-Order-Section For Stability
+%figure(1)
+%freqz(SOS, 4096, Fs) % Check Filter Performance
+%s_filtered = filtfilt(SOS, G, signal);
+
+fcutlow = 3000;
+fcuthigh = 3500;
+Wp = [fcutlow fcuthigh]/Fn; % Passband Frequency (Normalised)
+Ws = [fcutlow*0.95 fcuthigh/0.95]/Fn; % Stopband Frequency (Normalised)
+Rp = 30; % Passband Ripple (dB)
+Rs = 30; % Stopband Ripple (dB)
+%[n,Wn] = buttord(Wp,Ws,Rp,Rs); % Filter Order
+n=2;
+Wn=3200/Fn;
+[z,p,k] = butter(n, Wn, 'low'); % Filter Design
+[sosbp,gbp] = zp2sos(z, p, k); % Convert To Second-Order-Section For Stability
+%freqz(sosbp, 2^20, Fs) % Filter Bode Plot
+signal_filtered = filtfilt(sosbp, gbp, signal); % Filter Signal
+figure(2)
+subplot(2,1,1)
+plot(signal)
+subplot(2,1,2)
+plot(signal_filtered);
+
+
+
+% output
+audiowrite('restored.flac', signal_filtered, Fs);

TP3_8/trash/DiogoEliseu_TP3_8_cheb.m

@@ -0,0 +1,45 @@
+%% Inicialização do ambiente
+clear; close all; clc
+
+%% Exercício 8
+[y, Fs] = audioread("Canto1.mp3"); % Signal and Sampling Frequency
+s = y(:,1);
+L = length(s);
+
+Fn = Fs/2;                                              % Nyquist Frequency
+Wp = [150  5800]/Fn;                                    % Normalised Passband
+Ws = [ 50  6100]/Fn;                                    % Normalised Stopband
+Rp =  1;                                                % Passband Ripple (dB)
+Rs = 30;                                                % Stopband Ripple (dB)
+[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs);                         % Chebyshev Type II Order
+[b,a] = cheby2(n,Rs,Ws);                                % IIR Filter Coefficients
+[SOS,G] = tf2sos(b,a);                                  % Convert To Second-Order-Section For Stability
+figure(1)
+freqz(SOS, 4096, Fs)                                    % Check Filter Performance
+
+s_filtered = filtfilt(SOS, G, s);
+
+fcuts = [680 690 710 720  1190 1205  1210 1220  6000  6100];    % Frequency Vector
+mags =   [1 0 1 0 1 0];                                         % Magnitude (Defines Passbands & Stopbands)
+devs = [0.05  0.01  0.05  0.01 0.05  0.01];                     % Allowable Deviations
+[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,Fs);
+n = n + rem(n,2);
+hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'scale');
+figure(2)
+freqz(hh,1, 4096, Fs)                                           % Filter Bode Plot
+filt_sig = filtfilt(hh, 1, s);                                  % Filter Signal
+FTS = fft(s)/L;                                                 % FFT Of Original Signal
+FTFS = fft(filt_sig)/L;                                         % FFT Of Filtered Signal
+Fv = linspace(0, 1, fix(L/2)+1)*Fn;                             % FFT Frequency Vector
+Iv = 1:length(Fv);                                              % Index Vector
+[pks,Frs] = findpeaks(abs(FTS(Iv))*2, Fv, 'MinPeakHeight',0.1); % Find Pure Tone Frequencies
+figure(3)
+semilogy(Fv, (abs(FTS(Iv)))*2)
+hold on
+plot(Fv, (abs(FTFS(Iv)))*2)
+hold off
+grid
+axis([0  2500    ylim])
+
+% output
+audiowrite('restored.flac', s_filtered, Fs);

TP3_8/trash/DiogoEliseu_TP3_8_filterDesigner.m

@@ -0,0 +1,30 @@
+%% Inicialização do ambiente
+clear; close all; clc
+
+%% Exercício 8
+% Input
+[y, Fs] = audioread("Canto1.mp3"); % Signal and Sampling Frequency
+signal = y(:,1);
+Fn = Fs/2; % Nyquist Frequency (Hz)
+%Ts = 1/Fs;
+L = length(signal); % Signal Length
+
+% Find the noise
+FTsignal = fft(signal)/L;           % Fourier Transform
+Fv = linspace(0, 1, fix(L/2)+1)*Fn; % Frequency Vector
+Iv = 1:numel(Fv);                   % Index Vector
+figure(1)
+% x is Hz, y is amplitude
+plot(Fv, abs(FTsignal(Iv))*2)
+grid
+
+% Filter
+signal_filtered = filter(nd_designed_sexy,signal);
+figure(2)
+subplot(2,1,1)
+plot(signal)
+subplot(2,1,2)
+plot(signal_filtered);
+
+% output
+audiowrite('restored.flac', signal_filtered, Fs);

TP3_8/trash/designed_sexy.m

@@ -0,0 +1,23 @@
+function Hd = designed_sexy
+%DESIGNED_SEXY Returns a discrete-time filter object.
+
+% MATLAB Code
+% Generated by MATLAB(R) 9.8 and Signal Processing Toolbox 8.4.
+% Generated on: 28-Dec-2020 08:55:51
+
+% Butterworth Lowpass filter designed using FDESIGN.LOWPASS.
+
+% All frequency values are in Hz.
+Fs = 44100;  % Sampling Frequency
+
+Fpass = 3200;        % Passband Frequency
+Fstop = 6900;        % Stopband Frequency
+Apass = 1;           % Passband Ripple (dB)
+Astop = 80;          % Stopband Attenuation (dB)
+match = 'stopband';  % Band to match exactly
+
+% Construct an FDESIGN object and call its BUTTER method.
+h  = fdesign.lowpass(Fpass, Fstop, Apass, Astop, Fs);
+Hd = design(h, 'butter', 'MatchExactly', match);
+
+% [EOF]

TP3_8/trash/nd_designed_sexy.m

@@ -0,0 +1,27 @@
+function Hd = nd_designed_sexy
+%ND_DESIGNED_SEXY Returns a discrete-time filter object.
+
+% MATLAB Code
+% Generated by MATLAB(R) 9.8 and DSP System Toolbox 9.10.
+% Generated on: 28-Dec-2020 11:09:34
+
+% Butterworth Bandpass filter designed using FDESIGN.BANDPASS.
+
+% All frequency values are in Hz.
+Fs = 48000;  % Sampling Frequency
+
+Fstop1 = 730;         % First Stopband Frequency
+Fpass1 = 3000;        % First Passband Frequency
+Fpass2 = 5000;        % Second Passband Frequency
+Fstop2 = 6800;        % Second Stopband Frequency
+Astop1 = 50;          % First Stopband Attenuation (dB)
+Apass  = 0.0005;      % Passband Ripple (dB)
+Astop2 = 50;          % Second Stopband Attenuation (dB)
+match  = 'passband';  % 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]