From 262a0f2bbe2553f597b848fe1737b1a3753fb851 Mon Sep 17 00:00:00 2001 From: Priec Date: Sun, 7 Dec 2025 00:10:47 +0100 Subject: [PATCH] sem2 done --- semestralka2/analyza_c.m | 156 ++++++++++++++++ semestralka2/analyza_d.m | 124 +++++++++++++ semestralka2/analyza_experimental.m | 265 ++++++++++++++++++++++++++++ semestralka2/final.m | 198 +++++++++++++++++++++ semestralka2/final_cleaned.m | 252 ++++++++++++++++++++++++++ semestralka2/final_cleaned_backup.m | 251 ++++++++++++++++++++++++++ semestralka2/spd_sumu.m | 207 ++++++++++++++++++++++ 7 files changed, 1453 insertions(+) create mode 100644 semestralka2/analyza_c.m create mode 100644 semestralka2/analyza_d.m create mode 100644 semestralka2/analyza_experimental.m create mode 100644 semestralka2/final.m create mode 100644 semestralka2/final_cleaned.m create mode 100644 semestralka2/final_cleaned_backup.m create mode 100644 semestralka2/spd_sumu.m diff --git a/semestralka2/analyza_c.m b/semestralka2/analyza_c.m new file mode 100644 index 0000000..bb9b7b5 --- /dev/null +++ b/semestralka2/analyza_c.m @@ -0,0 +1,156 @@ +clc; +close all; +clear all; + +% === Načítanie vstupnej stopy === +[y, Fs] = audioread('sem2.wav'); +t = (0:length(y)-1) / Fs; + +figure; %1 +plot(t, y); +xlabel('Time'); +ylabel('Amplitude'); +title('Waveform sem2.wav'); +grid on; + +[y1, Fs1] = audioread('flac.wav'); +[y2, Fs2] = audioread('flac2.wav'); + +t1 = (0:length(y1)-1) / Fs1; +t2 = (0:length(y2)-1) / Fs2; + +figure; %2 +plot(t1, y1); +xlabel('Time [s]'); +ylabel('Amplitude'); +title('Waveform flac.wav'); +grid on; + +figure; %3 +plot(t2, y2); +xlabel('Time [s]'); +ylabel('Amplitude'); +title('Waveform flac2.wav'); +grid on; + +% === FFT a spektrá === +N1 = length(y1); +N2 = length(y2); + +X1 = fft(y1); +X2 = fft(y2); + +freq_shift1 = (-N1/2 : N1/2 - 1) * (Fs1 / N1); +freq_shift2 = (-N2/2 : N2/2 - 1) * (Fs2 / N2); + +figure; %4 +subplot(2,1,1); +plot(freq_shift1, real(fftshift(X1)), 'b'); +title('Reálna časť posunutého spektra flac.wav'); +xlabel('Frekvencia [Hz]'); +ylabel('Reálna hodnota'); +grid on; +subplot(2,1,2); +plot(freq_shift1, imag(fftshift(X1)), 'r'); +title('Imaginárna časť posunutého spektra flac.wav'); +xlabel('Frekvencia [Hz]'); +ylabel('Imaginárna hodnota'); +grid on; + +figure; %5 +subplot(2,1,1); +plot(freq_shift2, real(fftshift(X2)), 'b'); +title('Reálna časť posunutého spektra flac2.wav'); +xlabel('Frekvencia [Hz]'); +ylabel('Reálna hodnota'); +grid on; +subplot(2,1,2); +plot(freq_shift2, imag(fftshift(X2)), 'r'); +title('Imaginárna časť posunutého spektra flac2.wav'); +xlabel('Frekvencia [Hz]'); +ylabel('Imaginárna hodnota'); +grid on; + +% === Nájdeme 5 dominujúcich frekvencií z flac2.wav === +halfN2 = floor(N2/2); +realX2 = real(X2(1:halfN2)); +imagX2 = imag(X2(1:halfN2)); +freq2 = (0:halfN2-1)*(Fs2/N2); + +[pks2r, locs2r] = findpeaks(abs(realX2), 'NPeaks', 3, 'SortStr', 'descend'); +f0_2r = freq2(locs2r); +f0_2r(f0_2r < 2) = []; % odstráni DC + +[pks2i, locs2i] = findpeaks(abs(imagX2), 'NPeaks', 2, 'SortStr', 'descend'); +f0_2i = freq2(locs2i); +f0_2i(f0_2i < 2) = []; + +disp('Najvýraznejšie frekvenčné zložky flac2.wav - reálna časť:'); +disp(f0_2r'); +disp('Najvýraznejšie frekvenčné zložky flac2.wav - imaginárna časť:'); +disp(f0_2i'); + +% === IIR Notch Filter Design == +f_all = unique([f0_2r(:); f0_2i(:)]); +if length(f_all) > 5 + f_all = f_all(1:5); +end + +r = 0.995; % stabilita +sos_all = []; +gain_all = 1; + +for k = 1:length(f_all) + wo = f_all(k) / (Fs2 / 2); + bw = (1 - r); + [b_notch, a_notch] = iirnotch(wo, bw); + [sos, g] = tf2sos(b_notch, a_notch); + sos_all = [sos_all; sos]; %#ok + gain_all = gain_all * g; +end + +% === Frekvenčná charakteristika filtra === +figure; +freqz(sos_all, 1024, Fs2); +title('Frekvenčná charakteristika IIR Notch filtra (5 zložiek z flac2.wav)'); +xlabel('Frekvencia [Hz]'); +ylabel('Zosilnenie [dB]'); +grid on; + +% === Filtrovanie flac2.wav === +filtered_y2 = filtfilt(sos_all, gain_all, double(y2)); + +audiowrite('filtered_flac2.wav', filtered_y2, Fs2); + +% === FFT po filtrovaní flac2 === +X2_filt = fft(filtered_y2); +freq_shift2_filt = (-N2/2 : N2/2 - 1) * (Fs2 / N2); + +figure; +subplot(2,1,1); +plot(freq_shift2_filt, real(fftshift(X2_filt)), 'b'); +title('Reálna časť spektra po filtrovaní (filtered flac2.wav)'); +xlabel('Frekvencia [Hz]'); +ylabel('Reálna hodnota'); +grid on; + +subplot(2,1,2); +plot(freq_shift2_filt, imag(fftshift(X2_filt)), 'r'); +title('Imaginárna časť spektra po filtrovaní (filtered flac2.wav)'); +xlabel('Frekvencia [Hz]'); +ylabel('Imaginárna hodnota'); +grid on; + +% === Aplikácia filtra na pôvodný sem2.wav === +filtered_sem2 = filtfilt(sos_all, gain_all, double(y)); +audiowrite('filtered_sem2.wav', filtered_sem2, Fs); + +% === Vizualizácia filtrovaného sem2 === +figure; +plot(filtered_sem2); +title('Filtrovaný signál sem2.wav (filter z flac2.wav)'); +xlabel('Vzorky'); +ylabel('Amplitúda'); +grid on; + +disp('Filtrovanie hotové. Spočetné spektrá flac2 aj výsledný sem2 boli uložené.'); diff --git a/semestralka2/analyza_d.m b/semestralka2/analyza_d.m new file mode 100644 index 0000000..faf7e7f --- /dev/null +++ b/semestralka2/analyza_d.m @@ -0,0 +1,124 @@ +clc; +close all; +clear; + +%% === LOAD INPUTS === +[y, Fs] = audioread('sem2.wav'); +[y2, Fs2] = audioread('flac2.wav'); + +if size(y,2) > 1, y = mean(y,2); end +if size(y2,2) > 1, y2 = mean(y2,2); end + +t = (0:length(y)-1) / Fs; +t2 = (0:length(y2)-1) / Fs2; + +%% === VISUALIZE ORIGINAL WAVEFORMS === +figure; +subplot(2,1,1); +plot(t, y); +xlabel('Time [s]'); ylabel('Amplitude'); +title('sem2.wav – Original Signal'); +grid on; + +subplot(2,1,2); +plot(t2, y2, 'r'); +xlabel('Time [s]'); ylabel('Amplitude'); +title('flac2.wav – Noise Sample'); +grid on; + +%% === FREQUENCY ANALYSIS (FFT of noise) === +N2 = length(y2); +Y2 = fft(y2); +f2 = (0:N2/2-1) * Fs2 / N2; + +magY2 = abs(Y2(1:N2/2)); +figure; +plot(f2, 20*log10(magY2 + eps), 'r'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]'); +title('Noise Spectrum (flac2.wav)'); +xlim([0 Fs2/2]); grid on; + +%% === DETECT DOMINANT NOISE FREQUENCIES === +[pks, locs] = findpeaks(magY2, 'MinPeakHeight', max(magY2)*0.05, ... + 'MinPeakDistance', round(10*N2/Fs2)); +f_peaks = f2(locs); +f_peaks = f_peaks(f_peaks > 20 & f_peaks < Fs/2); +f_peaks = sort(f_peaks(1:min(5,length(f_peaks)))); +fprintf('Detected noise frequencies [Hz]: %.1f ', f_peaks); +fprintf('\n'); + +%% === DESIGN CASCADE OF NOTCH FILTERS (z-plane zeros/poles) === +r = 0.995; +sos_all = []; +g_all = 1; + +for k = 1:length(f_peaks) + w0 = 2*pi*f_peaks(k)/Fs; + z = [exp(1j*w0) exp(-1j*w0)]; + p = r*z; + b = real(poly(z)); + a = real(poly(p)); + b = b / (sum(b)/sum(a)); % normalize DC gain + [sos_i, g_i] = tf2sos(b, a); + sos_all = [sos_all; sos_i]; + g_all = g_all * g_i; +end + +%% === VISUALIZE FILTER FREQUENCY RESPONSE === +figure; +freqz(sos_all, 2048, Fs); +title('Cascade of Notch Filters (Pole-Zero Design)'); + +%% === APPLY FILTER TO sem2.wav === +y_notch = filtfilt(sos_all, g_all, y); +y_notch = y_notch / (max(abs(y_notch)) + eps) * 0.95; +audiowrite('sem2_notch_filtered.wav', y_notch, Fs); +fprintf('Saved: sem2_notch_filtered.wav\n'); + +%% === ADDITIONAL SPEECH BANDPASS (100–8000 Hz) === +[z_bp, p_bp, k_bp] = butter(4, [100 8000]/(Fs/2), 'bandpass'); +[sos_bp, g_bp] = zp2sos(z_bp, p_bp, k_bp); +y_final = filtfilt(sos_bp, g_bp, y_notch); +y_final = y_final / (max(abs(y_final)) + eps) * 0.95; +audiowrite('sem2_final.wav', y_final, Fs); +fprintf('Saved: sem2_final.wav (Notch + Bandpass)\n'); + +%% === TIME DOMAIN COMPARISON === +figure; +subplot(3,1,1); +plot(t, y, 'b'); title('Original sem2.wav'); +xlabel('Time [s]'); ylabel('Amplitude'); grid on; +subplot(3,1,2); +plot(t, y_notch, 'g'); title('After Notch Filtering'); +xlabel('Time [s]'); ylabel('Amplitude'); grid on; +subplot(3,1,3); +plot(t, y_final, 'r'); title('Final Output (Notch + Bandpass)'); +xlabel('Time [s]'); ylabel('Amplitude'); grid on; + +%% === FREQUENCY DOMAIN COMPARISON === +N = length(y); +f = (0:N/2-1)*Fs/N; +Y = fft(y); +Y_notch = fft(y_notch); +Y_final = fft(y_final); + +figure; +subplot(3,1,1); +plot(f, 20*log10(abs(Y(1:N/2))+eps), 'b'); +title('Original sem2.wav Spectrum'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]'); +xlim([0 2000]); grid on; +subplot(3,1,2); +plot(f, 20*log10(abs(Y_notch(1:N/2))+eps), 'g'); +title('After Notch Filtering Spectrum'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]'); +xlim([0 2000]); grid on; +subplot(3,1,3); +plot(f, 20*log10(abs(Y_final(1:N/2))+eps), 'r'); +title('After Notch + Bandpass Spectrum'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]'); +xlim([0 2000]); grid on; + +fprintf('\n Processing complete. Files saved:\n'); +fprintf(' sem2_notch_filtered.wav (Notch only)\n'); +fprintf(' sem2_final.wav (Notch + Bandpass)\n'); diff --git a/semestralka2/analyza_experimental.m b/semestralka2/analyza_experimental.m new file mode 100644 index 0000000..c1e6560 --- /dev/null +++ b/semestralka2/analyza_experimental.m @@ -0,0 +1,265 @@ +clc; +close all; +clear all; + +%% === STAGE 0: Load Audio Files === +[y, Fs] = audioread('sem2.wav'); +[y1, Fs1] = audioread('flac.wav'); +[y2, Fs2] = audioread('flac2.wav'); + +% Ensure mono +if size(y,2) > 1, y = mean(y,2); end +if size(y1,2) > 1, y1 = mean(y1,2); end +if size(y2,2) > 1, y2 = mean(y2,2); end + +fprintf('Loaded audio files:\n'); +fprintf(' sem2.wav: %d samples, Fs = %d Hz\n', length(y), Fs); +fprintf(' flac.wav: %d samples, Fs = %d Hz\n', length(y1), Fs1); +fprintf(' flac2.wav: %d samples, Fs = %d Hz\n', length(y2), Fs2); + +%% === STAGE 1: IIR Notch Filter === +fprintf('\n=== STAGE 1: Notch Filtering ===\n'); + +N2 = length(y2); +X2 = fft(y2); +halfN2 = floor(N2/2); +freq2 = (0:halfN2-1) * (Fs2/N2); + +magX2 = abs(X2(1:halfN2)); +[pks, locs] = findpeaks(magX2, 'NPeaks', 10, 'SortStr', 'descend', 'MinPeakDistance', 5); +f_peaks = freq2(locs); +f_peaks = f_peaks(f_peaks > 10); +if length(f_peaks) > 5 + f_peaks = f_peaks(1:5); +end + +fprintf('Detected tonal noise frequencies: '); +fprintf('%.1f Hz, ', f_peaks); +fprintf('\n'); + +r = 0.995; +sos_all = []; +gain_all = 1; + +for k = 1:length(f_peaks) + wo = f_peaks(k) / (Fs2/2); + if wo > 0 && wo < 1 + bw = 0.005; + [b_notch, a_notch] = iirnotch(wo, bw); + [sos, g] = tf2sos(b_notch, a_notch); + sos_all = [sos_all; sos]; + gain_all = gain_all * g; + end +end + +if ~isempty(sos_all) + y_notch = filtfilt(sos_all, gain_all, double(y)); + y2_notch = filtfilt(sos_all, gain_all, double(y2)); +else + y_notch = double(y); + y2_notch = double(y2); +end + +fprintf('Notch filtering complete.\n'); + +%% === STAGE 2: Spectral Subtraction === +fprintf('\n=== STAGE 2: Spectral Subtraction ===\n'); +params.frameLen = 2048; +params.hopLen = 512; +params.nfft = 2048; +params.alpha = 4.0; +params.beta = 0.002; +params.gamma = 1.0; +[y_specsub, ~] = spectral_subtraction(y_notch, y2_notch, Fs, params); +fprintf('Spectral subtraction complete.\n'); + +%% === STAGE 3: MMSE-LSA Enhancement === +fprintf('\n=== STAGE 3: MMSE-LSA Enhancement ===\n'); +params_mmse.frameLen = 2048; +params_mmse.hopLen = 512; +params_mmse.nfft = 2048; +params_mmse.alphaDD = 0.98; +params_mmse.minGain = 0.1; +params_mmse.noisePowSmooth = 0.7; +[y_mmse, ~] = mmse_lsa_denoise(y_specsub, y2_notch, Fs, params_mmse); +fprintf('MMSE-LSA enhancement complete.\n'); + +%% === STAGE 4: Post-processing === +fprintf('\n=== STAGE 4: Post-processing ===\n'); +[b_hp, a_hp] = butter(2, 80/(Fs/2), 'high'); +y_final = filtfilt(b_hp, a_hp, y_mmse); +y_final = y_final / (max(abs(y_final)) + eps) * 0.95; +fprintf('Post-processing complete.\n'); + +%% === Save outputs === +audiowrite('sem2_stage1_notch.wav', y_notch/max(abs(y_notch)+eps)*0.95, Fs); +audiowrite('sem2_stage2_specsub.wav', y_specsub/max(abs(y_specsub)+eps)*0.95, Fs); +audiowrite('sem2_stage3_mmse.wav', y_mmse/max(abs(y_mmse)+eps)*0.95, Fs); +audiowrite('sem2_final_cleaned.wav', y_final, Fs); + +y2_cleaned = spectral_subtraction(y2_notch, y2_notch, Fs2, params); +y2_cleaned = mmse_lsa_denoise(y2_cleaned, y2_notch, Fs2, params_mmse); +y2_cleaned = y2_cleaned / (max(abs(y2_cleaned)) + eps) * 0.95; +audiowrite('flac2_final_cleaned.wav', y2_cleaned, Fs2); +fprintf('\n✅ All files saved!\n'); + +%% === SIMPLE VISUALIZATION (same style as older system) === +t = (0:length(y)-1)/Fs; + +figure; +subplot(3,1,1); +plot(t, y); grid on; +title('Original sem2.wav'); +xlabel('Time [s]'); ylabel('Amplitude'); + +subplot(3,1,2); +plot(t, y_notch); grid on; +title('After Notch Filtering'); +xlabel('Time [s]'); ylabel('Amplitude'); + +subplot(3,1,3); +plot(t, y_final); grid on; +title('Final Cleaned sem2\_final\_cleaned.wav'); +xlabel('Time [s]'); ylabel('Amplitude'); + +figure; +subplot(3,1,1); +plot((0:length(y2)-1)/Fs2, y2); grid on; +title('Original flac2.wav (Noise Sample)'); +xlabel('Time [s]'); ylabel('Amplitude'); + +subplot(3,1,2); +plot((0:length(y2_notch)-1)/Fs2, y2_notch); grid on; +title('After Notch Filtering (Noise)'); +xlabel('Time [s]'); ylabel('Amplitude'); + +subplot(3,1,3); +plot((0:length(y2_cleaned)-1)/Fs2, y2_cleaned); grid on; +title('Cleaned flac2.wav (Noise Reference)'); +xlabel('Time [s]'); ylabel('Amplitude'); + +fprintf('\nProcessing summary:\n'); +fprintf(' - Stage 1 removed %d tonal frequencies\n', length(f_peaks)); +fprintf(' - Stage 2 applied spectral subtraction (α = %.1f)\n', params.alpha); +fprintf(' - Stage 3 applied MMSE-LSA (α_DD = %.2f)\n', params_mmse.alphaDD); +fprintf(' - Final result saved as sem2_final_cleaned.wav (best quality)\n'); + + +function [y_out, noise_psd] = spectral_subtraction(x, noise_ref, Fs, p) + x = x(:); + noise_ref = noise_ref(:); + frameLen = p.frameLen; hopLen = p.hopLen; nfft = p.nfft; + alpha = p.alpha; beta = p.beta; gamma = p.gamma; + win = hamming(frameLen, 'periodic'); + + [N_stft, ~, ~] = stft(noise_ref, Fs, 'Window', win, ... + 'OverlapLength', frameLen-hopLen, 'FFTLength', nfft); + noise_psd = mean(abs(N_stft).^2, 2); + noise_psd = movmean(noise_psd, 5); + + [X_stft, ~, ~] = stft(x, Fs, 'Window', win, ... + 'OverlapLength', frameLen-hopLen, 'FFTLength', nfft); + + [nFreq, nFrames] = size(X_stft); + Y_stft = zeros(size(X_stft)); + + for k = 1:nFrames + X_frame = X_stft(:, k); + X_mag = abs(X_frame); + X_phase = angle(X_frame); + X_pow = X_mag.^2; + snr_post = X_pow ./ (noise_psd + eps); + alpha_freq = alpha * (1 + 0.5 * (1 - min(snr_post, 4)/4)); + Y_pow = X_pow - alpha_freq .* noise_psd; + Y_pow = max(Y_pow, beta * X_pow); + Y_mag = sqrt(max(Y_pow, 0)); + Y_stft(:, k) = Y_mag .* exp(1j * X_phase); + end + + y_out = istft(Y_stft, Fs, 'Window', win, ... + 'OverlapLength', frameLen-hopLen, ... + 'FFTLength', nfft, 'ConjugateSymmetric', true); + y_out = y_out(1:min(length(y_out), length(x))); + if length(y_out) < length(x) + y_out = [y_out; zeros(length(x)-length(y_out),1)]; + end +end + + +function [y_out, G_all] = mmse_lsa_denoise(x, noise_ref, Fs, p) + x = x(:); noise_ref = noise_ref(:); + frameLen = p.frameLen; hopLen = p.hopLen; nfft = p.nfft; + alphaDD = p.alphaDD; Gmin = p.minGain; + win = hamming(frameLen, 'periodic'); + + [N_stft, ~, ~] = stft(noise_ref, Fs, 'Window', win, ... + 'OverlapLength', frameLen-hopLen, 'FFTLength', nfft); + lambda_d = mean(abs(N_stft).^2, 2); + lambda_d = movmean(lambda_d, 5); + + [X_stft, ~, ~] = stft(x, Fs, 'Window', win, ... + 'OverlapLength', frameLen-hopLen, 'FFTLength', nfft); + + [nFreq, nFrames] = size(X_stft); + Y_stft = zeros(size(X_stft)); + G_all = zeros(nFreq, nFrames); + xi_prev = ones(nFreq, 1); + + for k = 1:nFrames + X_frame = X_stft(:, k); + X_mag = abs(X_frame); + X_phase = angle(X_frame); + X_pow = X_mag.^2; + gamma_k = X_pow ./ (lambda_d + eps); + xi_ml = max(gamma_k - 1, 0); + xi = alphaDD * xi_prev + (1 - alphaDD) * xi_ml; + xi = max(xi, 1e-6); + v = xi .* gamma_k ./ (1 + xi); + G = zeros(nFreq,1); + for i = 1:nFreq + if v(i) < 0.01 + G(i) = xi(i)/(1+xi(i)); + else + G(i) = xi(i)/(1+xi(i)) * exp(0.5 * expint_approx(v(i))); + end + end + G = max(G, Gmin); + G = movmean(G, 3); + Y_mag = G .* X_mag; + Y_stft(:,k) = Y_mag .* exp(1j*X_phase); + G_all(:,k) = G; + xi_prev = xi; + end + + y_out = istft(Y_stft, Fs, 'Window', win, ... + 'OverlapLength', frameLen-hopLen, ... + 'FFTLength', nfft, 'ConjugateSymmetric', true); + y_out = y_out(1:min(length(y_out), length(x))); + if length(y_out) < length(x) + y_out = [y_out; zeros(length(x)-length(y_out),1)]; + end +end + + +function y = expint_approx(x) + if x < 1e-10 + y = -log(1e-10) - 0.5772; + elseif x < 1 + y = -log(x) - 0.5772156649; + term = x; + for n = 1:50 + y = y + term / n; + term = -term * x / (n + 1); + if abs(term) < 1e-10, break; end + end + else + y = exp(-x) / x; + term = 1; + for n = 1:10 + term = -term * n / x; + y = y + term * exp(-x) / x; + if abs(term) < 1e-10, break; end + end + end + y = -y; +end diff --git a/semestralka2/final.m b/semestralka2/final.m new file mode 100644 index 0000000..ba95515 --- /dev/null +++ b/semestralka2/final.m @@ -0,0 +1,198 @@ +clc; +close all; +clear; + +%% === LOAD AUDIO FILES === +[y1, Fs1] = audioread('flac.wav'); +[y2, Fs2] = audioread('flac2.wav'); + +if size(y1,2) > 1, y1 = mean(y1,2); end +if size(y2,2) > 1, y2 = mean(y2,2); end + +t1 = (0:length(y1)-1) / Fs1; +t2 = (0:length(y2)-1) / Fs2; + +fprintf('Loaded flac.wav: %d samples, Fs = %d Hz\n', length(y1), Fs1); +fprintf('Loaded flac2.wav: %d samples, Fs = %d Hz\n', length(y2), Fs2); + +%% === TIME-DOMAIN WAVEFORMS === +figure; +subplot(2,1,1); +plot(t1, y1); +xlabel('Time [s]'); +ylabel('Amplitude'); +title('flac.wav – Time Domain'); +grid on; + +subplot(2,1,2); +plot(t2, y2, 'r'); +xlabel('Time [s]'); +ylabel('Amplitude'); +title('flac2.wav – Time Domain'); +grid on; + +%% === FREQUENCY-DOMAIN ANALYSIS (FFT) === +N1 = length(y1); +N2 = length(y2); + +X1 = fft(y1); +X2 = fft(y2); + +freq_shift1 = (-N1/2 : N1/2 - 1) * (Fs1 / N1); +freq_shift2 = (-N2/2 : N2/2 - 1) * (Fs2 / N2); + +% --- Real / Imag components --- +figure; +subplot(2,1,1); +plot(freq_shift2, real(fftshift(X2))); +xlabel('Frequency [Hz]'); +ylabel('Real Part'); +title('flac2.wav – Real Part of FFT'); +grid on; +subplot(2,1,2); +plot(freq_shift2, imag(fftshift(X2)), 'r'); +xlabel('Frequency [Hz]'); +ylabel('Imag Part'); +title('flac2.wav – Imaginary Part of FFT'); +grid on; + +% --- Magnitude Spectrum --- +figure; +plot(freq_shift2, 20*log10(abs(fftshift(X2)) + eps), 'y'); +xlabel('Frequency [Hz]'); +ylabel('Magnitude [dB]'); +title('flac2.wav – Magnitude Spectrum'); +xlim([0 2000]); +grid on; + +%% === PEAK DETECTION TO FIND TONAL NOISE FREQUENCIES === +halfN2 = floor(N2/2); +f2 = (0:halfN2-1)*(Fs2/N2); +mag2 = abs(X2(1:halfN2)); + +[pks, locs] = findpeaks(mag2, 'MinPeakHeight', max(mag2)*0.1, ... + 'MinPeakDistance', round(10*N2/Fs2)); + +f_tonal = f2(locs); +f_tonal = f_tonal(f_tonal > 30 & f_tonal < Fs2/2); +f_tonal = sort(f_tonal(1:min(5, numel(f_tonal)))); + +fprintf('\nDetected tonal noise frequencies (from flac2.wav):\n'); +fprintf(' %.1f Hz\n', f_tonal); + +%% === NOTCH FILTER DESIGN (Z-domain zeros/poles) === +r = 0.995; % pole radius +sos_all = []; +g_all = 1; + +for i = 1:length(f_tonal) + f0 = f_tonal(i); + w0 = 2*pi*f0/Fs1; % use Fs1 so it matches flac.wav + z = [exp(1j*w0) exp(-1j*w0)]; % zeros at ±ω0 + p = r*z; % poles slightly inside circle + b = real(poly(z)); + a = real(poly(p)); + b = b / (sum(b)/sum(a)); % normalize gain + [sos_i, g_i] = tf2sos(b, a); + sos_all = [sos_all; sos_i]; + g_all = g_all * g_i; +end + +%% === VISUALIZE NOTCH FILTER RESPONSE === +figure; +freqz(sos_all, 2048, Fs1); +title('Notch Filter Cascade – Based on flac2 Noise Frequencies'); + +%% === APPLY NOTCH FILTERS TO SIGNALS === +y1_filtered = filtfilt(sos_all, g_all, y1); % Apply to flac.wav +y2_filtered = filtfilt(sos_all, g_all, y2); % Apply to flac2.wav +y1_filtered = y1_filtered / max(abs(y1_filtered)); +y2_filtered = y2_filtered / max(abs(y2_filtered)); + +audiowrite('filtered_flac.wav', y1_filtered, Fs1); +audiowrite('filtered_flac2.wav', y2_filtered, Fs2); + +fprintf('\nSaved:\n'); +fprintf(' filtered_flac.wav (filtered signal)\n'); +fprintf(' filtered_flac2.wav (filtered noise sample)\n'); + +%% === COMPARE ORIGINAL VS FILTERED SPECTRA === +N_plot = length(y1); +f_plot = (0:N_plot/2-1) * Fs1 / N_plot; +X1_orig = fft(y1); +X1_filt = fft(y1_filtered); + +figure; +subplot(2,1,1); +plot(f_plot, 20*log10(abs(X1_orig(1:N_plot/2))+eps)); +xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]'); +title('flac.wav – Original Spectrum'); +xlim([0 2000]); grid on; + +subplot(2,1,2); +plot(f_plot, 20*log10(abs(X1_filt(1:N_plot/2))+eps), 'r'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]'); +title('flac.wav – After Notch Filtering'); +xlim([0 2000]); grid on; + +%% === APPLY THE SAME FILTER TO sem2.wav === +fprintf('\n=== APPLYING DESIGNED FILTER TO sem2.wav ===\n'); + +[y_sem2, Fs_sem2] = audioread('sem2.wav'); +if size(y_sem2,2) > 1, y_sem2 = mean(y_sem2,2); end + +if Fs_sem2 ~= Fs1 + warning('Sample rates differ – resampling sem2.wav to match notch filter design'); + y_sem2 = resample(y_sem2, Fs1, Fs_sem2); + Fs_sem2 = Fs1; +end + +% --- Step 1: Apply multi‑notch filter from flac2 noise analysis --- +y_sem2_notch = filtfilt(sos_all, g_all, y_sem2); + +% --- Step 2 (optional): add band‑pass filter for speech 100‑8000 Hz --- +[z_bp, p_bp, k_bp] = butter(4, [100 8000]/(Fs_sem2/2), 'bandpass'); +[sos_bp, g_bp] = zp2sos(z_bp, p_bp, k_bp); +y_sem2_final = filtfilt(sos_bp, g_bp, y_sem2_notch); + +% --- Normalize both versions --- +y_sem2_notch = y_sem2_notch / (max(abs(y_sem2_notch)) + eps) * 0.95; +y_sem2_final = y_sem2_final / (max(abs(y_sem2_final)) + eps) * 0.95; + +% --- Save results --- +audiowrite('sem2_notch_filtered.wav', y_sem2_notch, Fs_sem2); +audiowrite('sem2_final.wav', y_sem2_final, Fs_sem2); +fprintf('Saved:\n'); +fprintf(' sem2_notch_filtered.wav (notch filters only)\n'); +fprintf(' sem2_final.wav (notch + speech band‑pass)\n'); + +%% === VISUAL CHECK – BEFORE vs AFTER === +t_sem2 = (0:length(y_sem2)-1)/Fs_sem2; +N_sem2 = length(y_sem2); +f_sem2 = (0:N_sem2/2-1)*Fs_sem2/N_sem2; +Y_orig = fft(y_sem2); +Y_filt = fft(y_sem2_final); + +figure; +subplot(2,1,1); +plot(t_sem2, y_sem2); +xlabel('Time [s]'); ylabel('Amplitude'); +title('sem2.wav – Original Time Domain'); grid on; + +subplot(2,1,2); +plot(t_sem2, y_sem2_final, 'r'); +xlabel('Time [s]'); ylabel('Amplitude'); +title('sem2\_final.wav – After Notch + Band‑Pass'); grid on; + +figure; +subplot(2,1,1); +plot(f_sem2, 20*log10(abs(Y_orig(1:N_sem2/2))+eps)); +xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]'); +title('sem2.wav – Original Spectrum'); +xlim([0 2000]); grid on; + +subplot(2,1,2); +plot(f_sem2, 20*log10(abs(Y_filt(1:N_sem2/2))+eps), 'r'); +xlabel('Frequency [Hz]'); ylabel('Magnitude [dB]'); +title('sem2\_final.wav – After Filtering'); +xlim([0 2000]); grid on; diff --git a/semestralka2/final_cleaned.m b/semestralka2/final_cleaned.m new file mode 100644 index 0000000..94933ed --- /dev/null +++ b/semestralka2/final_cleaned.m @@ -0,0 +1,252 @@ +clc; +close all; +clear; + +% === Global figure styling (light mode defaults) === +set(0, 'DefaultFigureColor', 'w'); % white figure background +set(0, 'DefaultAxesColor', 'w'); % white plotting area +set(0, 'DefaultAxesXColor', 'k'); % black x-axis lines and text +set(0, 'DefaultAxesYColor', 'k'); % black y-axis lines and text +set(0, 'DefaultTextColor', 'k'); % black text for titles, labels +set(0, 'DefaultLineColor', 'b'); % blue lines by default +clear saveFigure + +% === Automatic figure export counter === +fig_counter = 1; +save_fig = @saveFigure; + +function saveFigure(name) + persistent counter + if isempty(counter) + counter = 1; + end + + % Force painters renderer and white backgrounds + set(gcf, 'Renderer', 'painters'); + set(gcf, 'Color', [1 1 1]); + axesHandles = findall(gcf, 'Type', 'axes'); + for ax = axesHandles' + set(ax, 'Color', [1 1 1]); + end + + % Export clean EPS with no transparency problems + print(gcf, '-depsc2', '-r300', '-painters', sprintf('%d_%s.eps', counter, name)); + + counter = counter + 1; +end + +[y1, Fs1] = audioread('flac.wav'); +[y2, Fs2] = audioread('flac2.wav'); +[y_sem2, Fs_sem2] = audioread('sem2.wav'); + +% Časové priebehy +t1 = (0:length(y1)-1)/Fs1; +t2 = (0:length(y2)-1)/Fs2; +t_sem2 = (0:length(y_sem2)-1)/Fs_sem2; + +figPos = get(0, 'DefaultFigurePosition'); +figPos(4) = figPos(4) / 2; % polovičná výška (4 = height element) +figure('Position', figPos); % style defaults +plot(t_sem2, y_sem2); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('sem2.wav – pôvodný (RAW) signál'); +grid on; +save_fig('orig_sem2'); fig_counter = fig_counter + 1; + +figure; +subplot(2,1,1); +plot(t1, y1); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('flac.wav'); +grid on; + +subplot(2,1,2); +plot(t2, y2, 'r'); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('flac2.wav'); +grid on; +save_fig('time'); fig_counter = fig_counter + 1; + +% Fourierova transformácia +N1 = length(y1); +N2 = length(y2); +X1 = fft(y1); +X2 = fft(y2); + +freq_shift1 = (-N1/2 : N1/2 - 1) * (Fs1 / N1); +freq_shift2 = (-N2/2 : N2/2 - 1) * (Fs2 / N2); + +% Real a imag časť spektra +figure; +subplot(2,1,1); +plot(freq_shift2, real(fftshift(X2))); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac2.wav – Reálna časť spektra'); +xlim([-550 550]); +grid on; + +subplot(2,1,2); +plot(freq_shift2, imag(fftshift(X2)), 'r'); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac2.wav – Imaginárna časť spektra'); +xlim([-550 550]); +grid on; +save_fig('spectrum'); fig_counter = fig_counter + 1; + +% Real a imag časť spektra +figure; +subplot(2,1,1); +plot(freq_shift1, real(fftshift(X1))); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac.wav – Reálna časť spektra'); +xlim([-3500 3500]); +grid on; + +subplot(2,1,2); +plot(freq_shift1, imag(fftshift(X1)), 'r'); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac.wav – Imaginárna časť spektra'); +xlim([-3500 3500]); +grid on; +save_fig('spectrum2'); fig_counter = fig_counter + 1; + +% Detekcia dominantných frekvencií +halfN2 = floor(N2 / 2); +freq = (0:halfN2-1) * (Fs2 / N2); +magX = abs(X2(1:halfN2)); + +[pks, locs] = findpeaks(magX, 'NPeaks', 3, 'SortStr', 'descend'); +f_tonal = sort(freq(locs)); + +fprintf('\nDetekovane frekvencie z (flac2.wav):\n'); +fprintf(' %.1f Hz\n', f_tonal); + +% SOS IIR filter +r = 0.99; +Fs = Fs1; +sos = zeros(length(f_tonal), 6); + +for k = 1:length(f_tonal) + theta = 2 * pi * (f_tonal(k) / Fs); + b = [1, -2*cos(theta), 1]; + a = [1, -2*r*cos(theta), r^2]; + sos(k, :) = [b, a]; +end + +% Nuly a póly +figure; +[z, p, k_gain] = sos2zp(sos); +zplane(z, p); +title('Nuly a póly notch filtrov'); +xlabel('Reálna časť'); +ylabel('Imaginárna časť'); +save_fig('pz_notch'); fig_counter = fig_counter + 1; + +% Frekvenčná charakteristika +figure; +[h, w] = freqz(sos, 4096, Fs); +freqz(sos, 4096, Fs); +title('Frekvenčná charakteristika notch filtra'); +subplot(2,1,1); xlim([0 0.5]); +subplot(2,1,2); xlim([0 0.5]); +save_fig('freqz_notch'); fig_counter = fig_counter + 1; + +% Filtrovanie flac.wav a flac2.wav +y1_filtered = sosfilt(sos, double(y1)); +y2_filtered = sosfilt(sos, double(y2)); + +figure; +subplot(2,1,1); +plot(y2); +title('Pôvodný signál flac2.wav'); +xlabel('Vzorky'); +ylabel('Amplitúda'); +ax1 = gca; +grid on; + +subplot(2,1,2); +plot(y2_filtered, 'r'); +title('Filtrovaný signál flac2.wav'); +xlabel('Vzorky'); +ylabel('Amplitúda'); +ax2 = gca; +grid on; + +% spolocna y os pre oba grafy +ylim_common = [min([ax1.YLim, ax2.YLim]), max([ax1.YLim, ax2.YLim])]; +ax1.YLim = ylim_common; +ax2.YLim = ylim_common; +save_fig('orig_vs_filtered_flac2'); fig_counter = fig_counter + 1; + +% Aplikácia filtrov na sem2.wav +y_sem2_filtered = sosfilt(sos, double(y_sem2)); + +figure; +subplot(2,1,1); +plot(t_sem2, y_sem2); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('sem2.wav – pôvodný signál'); +grid on; + +subplot(2,1,2); +plot(t_sem2, y_sem2_filtered, 'r'); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('sem2.wav – filtrovaný signál'); +grid on; + +ylim_common = [min([min(y_sem2), min(y_sem2_filtered)]), ... + max([max(y_sem2), max(y_sem2_filtered)])]; +subplot(2,1,1); ylim(ylim_common); +subplot(2,1,2); ylim(ylim_common); +save_fig('orig_vs_filtered_sem2'); fig_counter = fig_counter + 1; + +audiowrite('filtered_sem2.wav', y_sem2_filtered, Fs_sem2); + +% Butterworth filter pre rozsah ľudskej reči +[z_bp, p_bp, k_bp] = butter(4, [100 8000]/(Fs_sem2/2), 'bandpass'); +[sos_bp, g_bp] = zp2sos(z_bp, p_bp, k_bp); +y_sem2_final = filtfilt(sos_bp, g_bp, y_sem2_filtered); + +audiowrite('filtered_sem2_voiceband.wav', ... + y_sem2_final / max(abs(y_sem2_final)), Fs_sem2); + +t_sem2 = (0:length(y_sem2_filtered)-1) / Fs_sem2; + +figure; +subplot(2,1,1); +plot(t_sem2, y_sem2_filtered); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('sem2.wav – notch filtrácia'); +grid on; + +subplot(2,1,2); +plot(t_sem2, y_sem2_final, 'r'); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('sem2.wav – notch + band-pass filtrácia'); +grid on; + +ylim_common = [min([min(y_sem2_filtered), min(y_sem2_final)]), ... + max([max(y_sem2_filtered), max(y_sem2_final)])]; +subplot(2,1,1); ylim(ylim_common); +subplot(2,1,2); ylim(ylim_common); +save_fig('notch_vs_bandpass_sem2'); fig_counter = fig_counter + 1; + +% Percentuálny rozdiel band-pass filtrácie +diff_signal = y_sem2_final - y_sem2_filtered; +power_original = sum(y_sem2_filtered.^2); +power_diff = sum(diff_signal.^2); +percent_difference = (power_diff / power_original) * 100; +fprintf('\nRozdiel medzi notch a notch+band‑pass signálom: %.2f %%\n', percent_difference); + + diff --git a/semestralka2/final_cleaned_backup.m b/semestralka2/final_cleaned_backup.m new file mode 100644 index 0000000..40253f4 --- /dev/null +++ b/semestralka2/final_cleaned_backup.m @@ -0,0 +1,251 @@ +clc; +close all; +clear; + +% === Global figure styling (light mode defaults) === +set(0, 'DefaultFigureColor', 'w'); % white figure background +set(0, 'DefaultAxesColor', 'w'); % white plotting area +set(0, 'DefaultAxesXColor', 'k'); % black x-axis lines and text +set(0, 'DefaultAxesYColor', 'k'); % black y-axis lines and text +set(0, 'DefaultTextColor', 'k'); % black text for titles, labels +set(0, 'DefaultLineColor', 'b'); % blue lines by default +clear saveFigure + +% === Automatic figure export counter === +fig_counter = 1; +save_fig = @saveFigure; + +function saveFigure(name) + persistent counter + if isempty(counter) + counter = 1; + end + + % Force painters renderer and white backgrounds + set(gcf, 'Renderer', 'painters'); + set(gcf, 'Color', [1 1 1]); + axesHandles = findall(gcf, 'Type', 'axes'); + for ax = axesHandles' + set(ax, 'Color', [1 1 1]); + end + + % Export clean EPS with no transparency problems + print(gcf, '-depsc2', '-r300', '-painters', sprintf('%d_%s.eps', counter, name)); + + counter = counter + 1; +end + +[y1, Fs1] = audioread('flac.wav'); +[y2, Fs2] = audioread('flac2.wav'); + +% Časové priebehy +t1 = (0:length(y1)-1)/Fs1; +t2 = (0:length(y2)-1)/Fs2; + +figure; +subplot(2,1,1); +plot(t1, y1); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('flac.wav'); +grid on; + +subplot(2,1,2); +plot(t2, y2, 'r'); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('flac2.wav'); +grid on; +save_fig('time'); fig_counter = fig_counter + 1; + +% Fourierova transformácia +N1 = length(y1); +N2 = length(y2); +X1 = fft(y1); +X2 = fft(y2); + +freq_shift1 = (-N1/2 : N1/2 - 1) * (Fs1 / N1); +freq_shift2 = (-N2/2 : N2/2 - 1) * (Fs2 / N2); + +% Real a imag časť spektra +figure; +subplot(2,1,1); +plot(freq_shift2, real(fftshift(X2))); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac2.wav – Reálna časť spektra'); +xlim([-550 550]); +grid on; + +subplot(2,1,2); +plot(freq_shift2, imag(fftshift(X2)), 'r'); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac2.wav – Imaginárna časť spektra'); +xlim([-550 550]); +grid on; +save_fig('spectrum'); fig_counter = fig_counter + 1; + +% Real a imag časť spektra +figure; +subplot(2,1,1); +plot(freq_shift1, real(fftshift(X1))); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac.wav – Reálna časť spektra'); +xlim([-3500 3500]); +grid on; + +subplot(2,1,2); +plot(freq_shift1, imag(fftshift(X1)), 'r'); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac.wav – Imaginárna časť spektra'); +xlim([-3500 3500]); +grid on; +save_fig('spectrum2'); fig_counter = fig_counter + 1; + +% Detekcia dominantných frekvencií +halfN2 = floor(N2 / 2); +freq = (0:halfN2-1) * (Fs2 / N2); +magX = abs(X2(1:halfN2)); + +[pks, locs] = findpeaks(magX, 'NPeaks', 3, 'SortStr', 'descend'); +f_tonal = sort(freq(locs)); + +fprintf('\nDetekovane frekvencie z (flac2.wav):\n'); +fprintf(' %.1f Hz\n', f_tonal); + +% SOS IIR filter +r = 0.99; +Fs = Fs1; +sos = zeros(length(f_tonal), 6); + +for k = 1:length(f_tonal) + theta = 2 * pi * (f_tonal(k) / Fs); + b = [1, -2*cos(theta), 1]; + a = [1, -2*r*cos(theta), r^2]; + sos(k, :) = [b, a]; +end + +% Nuly a póly +figure; +[z, p, k_gain] = sos2zp(sos); +zplane(z, p); +title('Nuly a póly notch filtrov'); +xlabel('Reálna časť'); +ylabel('Imaginárna časť'); +save_fig('pz_notch'); fig_counter = fig_counter + 1; + +% Frekvenčná charakteristika +figure; +[h, w] = freqz(sos, 4096, Fs); +freqz(sos, 4096, Fs); +title('Frekvenčná charakteristika notch filtra'); +subplot(2,1,1); xlim([0 0.5]); +subplot(2,1,2); xlim([0 0.5]); +save_fig('freqz_notch'); fig_counter = fig_counter + 1; + +% Filtrovanie flac.wav a flac2.wav +y1_filtered = sosfilt(sos, double(y1)); +y2_filtered = sosfilt(sos, double(y2)); + +figure; +subplot(2,1,1); +plot(y2); +title('Pôvodný signál flac2.wav'); +xlabel('Vzorky'); +ylabel('Amplitúda'); +ax1 = gca; +grid on; + +subplot(2,1,2); +plot(y2_filtered, 'r'); +title('Filtrovaný signál flac2.wav'); +xlabel('Vzorky'); +ylabel('Amplitúda'); +ax2 = gca; +grid on; + +ylim_common = [min([ax1.YLim, ax2.YLim]), max([ax1.YLim, ax2.YLim])]; +ax1.YLim = ylim_common; +ax2.YLim = ylim_common; +save_fig('orig_vs_filtered_flac2'); fig_counter = fig_counter + 1; + +% Aplikácia filtrov na sem2.wav +[y_sem2, Fs_sem2] = audioread('sem2.wav'); +y_sem2_filtered = sosfilt(sos, double(y_sem2)); + +t_sem2 = (0:length(y_sem2)-1) / Fs_sem2; + +figure; +subplot(2,1,1); +plot(t_sem2, y_sem2); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('sem2.wav – pôvodný signál'); +grid on; + +subplot(2,1,2); +plot(t_sem2, y_sem2_filtered, 'r'); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('sem2.wav – filtrovaný signál'); +grid on; + +ylim_common = [min([min(y_sem2), min(y_sem2_filtered)]), ... + max([max(y_sem2), max(y_sem2_filtered)])]; +subplot(2,1,1); ylim(ylim_common); +subplot(2,1,2); ylim(ylim_common); +save_fig('orig_vs_filtered_sem2'); fig_counter = fig_counter + 1; + +audiowrite('filtered_sem2.wav', y_sem2_filtered, Fs_sem2); + +% Butterworth filter pre rozsah ľudskej reči +[z_bp, p_bp, k_bp] = butter(4, [100 8000]/(Fs_sem2/2), 'bandpass'); +[sos_bp, g_bp] = zp2sos(z_bp, p_bp, k_bp); +y_sem2_final = filtfilt(sos_bp, g_bp, y_sem2_filtered); + +audiowrite('filtered_sem2_voiceband.wav', ... + y_sem2_final / max(abs(y_sem2_final)), Fs_sem2); + +t_sem2 = (0:length(y_sem2_filtered)-1) / Fs_sem2; + +figure; +subplot(2,1,1); +plot(t_sem2, y_sem2_filtered); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('sem2.wav – notch filtrácia'); +grid on; + +subplot(2,1,2); +plot(t_sem2, y_sem2_final, 'r'); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('sem2.wav – notch + band-pass filtrácia'); +grid on; + +ylim_common = [min([min(y_sem2_filtered), min(y_sem2_final)]), ... + max([max(y_sem2_filtered), max(y_sem2_final)])]; +subplot(2,1,1); ylim(ylim_common); +subplot(2,1,2); ylim(ylim_common); +save_fig('notch_vs_bandpass_sem2'); fig_counter = fig_counter + 1; + +% Percentuálny rozdiel band-pass filtrácie +diff_signal = y_sem2_final - y_sem2_filtered; +power_original = sum(y_sem2_filtered.^2); +power_diff = sum(diff_signal.^2); +percent_difference = (power_diff / power_original) * 100; +fprintf('\nRozdiel medzi notch a notch+band‑pass signálom: %.2f %%\n', percent_difference); + +figPos = get(0, 'DefaultFigurePosition'); +figPos(4) = figPos(4) / 2; % polovičná výška (4 = height element) +figure('Position', figPos); % použije globálne defaulty štýlu +plot(t_sem2, y_sem2); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('sem2.wav – pôvodný (RAW) signál'); +grid on; +save_fig('orig_sem2'); fig_counter = fig_counter + 1; + diff --git a/semestralka2/spd_sumu.m b/semestralka2/spd_sumu.m new file mode 100644 index 0000000..efa499c --- /dev/null +++ b/semestralka2/spd_sumu.m @@ -0,0 +1,207 @@ +clc; +close all; +clear; + +% === Global figure styling (light mode defaults) === +set(0, 'DefaultFigureColor', 'w'); % white figure background +set(0, 'DefaultAxesColor', 'w'); % white plotting area +set(0, 'DefaultAxesXColor', 'k'); % black x-axis lines and text +set(0, 'DefaultAxesYColor', 'k'); % black y-axis lines and text +set(0, 'DefaultTextColor', 'k'); % black text for titles, labels +set(0, 'DefaultLineColor', 'b'); % blue lines by default +clear saveFigure + +% === Automatic figure export counter === +fig_counter = 1; +save_fig = @saveFigure; + +function saveFigure(name) + persistent counter + if isempty(counter) + counter = 1; + end + + % Force painters renderer and white backgrounds + set(gcf, 'Renderer', 'painters'); + set(gcf, 'Color', [1 1 1]); + axesHandles = findall(gcf, 'Type', 'axes'); + for ax = axesHandles' + set(ax, 'Color', [1 1 1]); + end + + % Export clean EPS with no transparency problems + print(gcf, '-depsc2', '-r300', '-painters', sprintf('%d_%s.eps', counter, name)); + + counter = counter + 1; +end + +[y1, Fs1] = audioread('flac.wav'); +[y2, Fs2] = audioread('flac2.wav'); + +% Časové priebehy +t1 = (0:length(y1)-1)/Fs1; +t2 = (0:length(y2)-1)/Fs2; + +figure; +subplot(2,1,1); +plot(t1, y1); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('flac.wav'); +grid on; + +subplot(2,1,2); +plot(t2, y2, 'r'); +xlabel('Čas [s]'); +ylabel('Amplitúda'); +title('flac2.wav'); +grid on; +save_fig('time'); fig_counter = fig_counter + 1; + +% Fourierova transformácia +N1 = length(y1); +N2 = length(y2); +X1 = fft(y1); +X2 = fft(y2); + +freq_shift1 = (-N1/2 : N1/2 - 1) * (Fs1 / N1); +freq_shift2 = (-N2/2 : N2/2 - 1) * (Fs2 / N2); + +% Real a imag časť spektra +figure; +subplot(2,1,1); +plot(freq_shift2, real(fftshift(X2))); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac2.wav – Reálna časť spektra'); +xlim([-550 550]); +grid on; + +subplot(2,1,2); +plot(freq_shift2, imag(fftshift(X2)), 'r'); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac2.wav – Imaginárna časť spektra'); +xlim([-550 550]); +grid on; +save_fig('spectrum'); fig_counter = fig_counter + 1; + +% Real a imag časť spektra +figure; +subplot(2,1,1); +plot(freq_shift1, real(fftshift(X1))); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac.wav – Reálna časť spektra'); +xlim([-3500 3500]); +grid on; + +subplot(2,1,2); +plot(freq_shift1, imag(fftshift(X1)), 'r'); +xlabel('Frekvencia [Hz]'); +ylabel('Amplitúda'); +title('flac.wav – Imaginárna časť spektra'); +xlim([-3500 3500]); +grid on; +save_fig('spectrum2'); fig_counter = fig_counter + 1; + +% Určenie typu šumu na základe spektrálneho sklonu +[PSD, f] = pwelch(y1, hamming(4096), [], [], Fs1, 'onesided'); % PSD +mask = f > 10 & f < Fs1/2; +logf = log10(f(mask)); +logP = log10(PSD(mask)); +p = polyfit(logf, logP, 1); +slope = p(1); +intercept = p(2); +fprintf('Spektrálny sklon (alfa) = %.2f\n', slope); + +% Vizualizácia +figure; +loglog(f, PSD); hold on; +loglog(f(mask), 10.^(polyval(p, logf)), 'r', 'LineWidth', 1.5); +xlabel('Frekvencia [Hz]'); +ylabel('Výkonová spektrálna hustota (PSD)'); +title(sprintf('Spektrálny sklon = %.2f', slope)); +lgd = legend('Meraná PSD', 'Lineárna aproximácia (log-log)', 'Location', 'best'); +set(lgd, 'Color', 'w', 'TextColor', 'k'); +grid on; +save_fig('spectral_slope'); fig_counter = fig_counter + 1; + +% === IIR Notch filtering flac2.wav and noise type analysis === + +fprintf('\n=== IIR Notch Filtering of flac2.wav ===\n'); + +% Detekcia dominantných frekvencií v flac2.wav pre notch filtráciu +halfN = floor(N2 / 2); +freq_pos = (0:halfN - 1) * (Fs2 / N2); +magX2 = abs(X2(1:halfN)); +[pks2, locs2] = findpeaks(magX2, 'NPeaks', 3, 'SortStr', 'descend'); +f_notch = sort(freq_pos(locs2)); + +fprintf('\nDetekované dominantné frekvencie vo flac2.wav:\n'); +fprintf(' %.1f Hz\n', f_notch); + +% === Konštrukcia a aplikácia IIR notch filtra === +r = 0.98; % ostrosť zárezu (0.98 až 0.995) +Fs = Fs2; +sos_notch = zeros(length(f_notch), 6); + +for k = 1:length(f_notch) + theta = 2 * pi * (f_notch(k) / Fs); + b = [1, -2*cos(theta), 1]; + a = [1, -2*r*cos(theta), r^2]; + sos_notch(k, :) = [b, a]; +end + +% Aplikácia filtra na flac2.wav +y2_notched = sosfilt(sos_notch, double(y2)); + +% Uloženie výsledného signálu +audiowrite('flac2_notched.wav', y2_notched / max(abs(y2_notched)), Fs2); + +% Zobrazenie pred/po filtrácii +figure; +subplot(2,1,1); +plot(y2); +xlabel('Vzorky'); ylabel('Amplitúda'); +title('Pôvodný signál flac2.wav'); +grid on; + +subplot(2,1,2); +plot(y2_notched, 'r'); +xlabel('Vzorky'); ylabel('Amplitúda'); +title('Filtrovaný signál – flac2.wav (IIR Notch)'); +grid on; +save_fig('flac2_notched'); fig_counter = fig_counter + 1; + +% === Analýza zvyškového šumu po IIR filtri === +[PSD_res, f_res] = pwelch(y2_notched, hamming(4096), [], [], Fs2, 'onesided'); +mask_res = f_res > 10 & f_res < Fs2/2; +logf_res = log10(f_res(mask_res)); +logP_res = log10(PSD_res(mask_res)); +p_res = polyfit(logf_res, logP_res, 1); +slope_res = p_res(1); + +fprintf('\nSpektrálny sklon po IIR filtrácii (alfa) = %.2f\n', slope_res); + +% Vizualizácia +figure; +loglog(f_res, PSD_res); hold on; +loglog(f_res(mask_res), 10.^(polyval(p_res, logf_res)), 'r', 'LineWidth', 1.5); +xlabel('Frekvencia [Hz]'); +ylabel('Výkonová spektrálna hustota (PSD)'); +title(sprintf('Zvyškový šum po IIR notch filtrácii (α = %.2f)', slope_res)); +legend('PSD po IIR filtri', 'Lineárna aproximácia', 'Location', 'best'); +grid on; +save_fig('residual_noise_flac2'); fig_counter = fig_counter + 1; + +% Interpretácia typu šumu podľa alfa +if slope_res > -0.5 + noise_type = 'biely šum'; +elseif slope_res > -1.5 + noise_type = 'ružový šum'; +else + noise_type = 'hnedý šum'; +end + +fprintf('Typ zostávajúceho šumu po notch filtrácii: %s\n', noise_type);