From e1af76be20dee4c2f6b6f7ab34032e2084cc7ebc Mon Sep 17 00:00:00 2001
From: DUBOC Marc <marc.duboc@imt-atlantique.net>
Date: Fri, 2 May 2025 08:46:40 +0000
Subject: [PATCH] Upload New File
---
PARTIE3.m | 189 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 189 insertions(+)
create mode 100644 PARTIE3.m
diff --git a/PARTIE3.m b/PARTIE3.m
new file mode 100644
index 0000000..888691c
--- /dev/null
+++ b/PARTIE3.m
@@ -0,0 +1,189 @@
+%% ========================================================================
+% PARTIE 3 — Effets audio numeriques (Q 3.1 → Q 3.19)
+% Poly « SAR – Traitements audio : partie signal » v1.0 (23 avril 2025)
+% ------------------------------------------------------------------------
+clear; close all; clc
+Fe = 44100; % frequence d echantillonnage
+
+%% ========================================================================
+%% Q 3.1 — Correlation croisee R yx = R xx * h
+% -------------------------------------------------------------------------
+dur = 0.05; f0 = 1000;
+t = (0:1/Fe:dur - 1/Fe).';
+x = sin(2*pi*f0*t) + 0.5*randn(size(t)); % excitation
+h = [1 zeros(1,199) 0.5]; % reponse impulsionnelle test
+y = conv(x,h,'same');
+
+lagMax = 400;
+Rxx = xcorr(x,lagMax,'biased');
+Ryx = xcorr(y,x,lagMax,'biased'); % variable corrigee
+convTest = conv(Rxx, h, 'same'); % theorie
+
+figure('Name','Q3-1 Correlations'); hold on
+plot(-lagMax:lagMax, Ryx,'b', 'DisplayName','R yx numerique');
+plot(-lagMax:lagMax, convTest,'r--','DisplayName','R xx * h');
+grid on; legend; xlabel('retard (ech)');
+title('Q 3.1 Verification corrélation');
+
+%% ========================================================================
+%% Q 3.2 — Estimation de h si R xx ≈ Dirac
+% -------------------------------------------------------------------------
+hEstime = Ryx(lagMax+1:end); % partie causale
+figure('Name','Q3-2 h estimee');
+stem(hEstime(1:300)); grid on
+title('Q 3.2 Estimation h');
+
+%% ========================================================================
+%% Q 3.3 — Choix du meilleur signal d excitation
+% -------------------------------------------------------------------------
+if exist('signalExcitation.mat','file')
+ load signalExcitation.mat % attend variables xe1, xe2
+ lag = 1024;
+ Rxe1 = max(abs(xcorr(xe1,lag,'biased')));
+ Rxe2 = max(abs(xcorr(xe2,lag,'biased')));
+ if Rxe1 > Rxe2, choix = 'xe1'; else, choix = 'xe2'; end
+ fprintf('Q 3.3 Signal retenu : %s\n', choix);
+else
+ warning('Fichier signalExcitation.mat absent — Q 3.3 ignoree');
+ choix = '';
+end
+
+%% ========================================================================
+%% Q 3.4 — Simulation de mesure de piece et estimation de h
+% -------------------------------------------------------------------------
+if exist('simulePiece','file') && ~isempty(choix)
+ [tensionMicro, excitation] = simulePiece(choix);
+ lag = 4000;
+ hPiece = xcorr(tensionMicro, excitation, lag,'biased');
+ hPiece = hPiece(lag+1:end);
+ figure('Name','Q3-4 h piece');
+ plot(hPiece); grid on; title('Q 3.4 h mesuree');
+else
+ warning('Fonction simulePiece absente ou signal non choisi — Q 3.4 ignoree');
+ hPiece = h; % secours : reponse test
+ excitation = x; % secours : signal test
+end
+
+%% ========================================================================
+%% Q 3.5 — Fonction effectReverb (convolution directe)
+% -------------------------------------------------------------------------
+yReverb = effectReverb(excitation, hPiece); %#ok<NASGU>
+
+%% ========================================================================
+%% Q 3.6 — Temps de calcul de effectReverb
+tic
+effectReverb(excitation, hPiece);
+tempsDirect = toc;
+fprintf('Q 3.6 Convolution temporelle : %.3f s\n', tempsDirect);
+
+%% ========================================================================
+%% Q 3.7 — Fonction effectReverbFFT + temps
+tic
+yReverbFFT = effectReverbFFT(excitation, hPiece);
+tempsFFT = toc;
+fprintf('Q 3.7 Convolution FFT : %.3f s\n', tempsFFT);
+
+%% ========================================================================
+%% Q 3.8 — Equivalence numerique des deux convolutions
+delta = yReverbFFT(1:length(excitation)) - yReverb(1:length(excitation));
+fprintf('Q 3.8 Energie difference : %.4e\n', norm(delta));
+
+%% ========================================================================
+%% Q 3.9 – Q 3.14 — Delay a boucle simple
+% -------------------------------------------------------------------------
+g = 0.9; tauSec = 0.25; tauEch = round(tauSec * Fe);
+Nimp = 8; % nombre d echos
+hDelayTheo = impulseDelay(g, tauEch, Nimp);
+[bDelay, aDelay] = coeffDelay(g, tauEch);
+deltaSig = [1 zeros(1, Nimp*tauEch)];
+hDelayNum = filter(bDelay, aDelay, deltaSig);
+
+% Q 3.12 : affichage des impulsions
+figure('Name','Q3-12 h delay'); hold on
+stem(hDelayNum,'b'); stem(hDelayTheo,'r.');
+legend('numerique','theorique'); title('Impulsion delay');
+
+% Q 3.13 / 3.14 : module theorique vs DFT
+Nfft = 4096;
+w = linspace(0, pi, Nfft);
+Hth = 1 ./ (1 + g * exp(-1j * w * tauEch));
+Hnum = fft(hDelayNum, Nfft);
+figure('Name','Q3-14 Module'); hold on
+plot(w/pi*Fe/2, 20*log10(abs(Hth)),'r','DisplayName','theorie');
+plot(w/pi*Fe/2, 20*log10(abs(Hnum(1:Nfft))),'b--','DisplayName','DFT');
+legend; grid on; xlabel('f (Hz)'); ylabel('dB');
+title('Module theorique vs numerique');
+
+%% ========================================================================
+%% Q 3.15 — Fonction effectDelay + test
+testDur = 2; t = (0:1/Fe:testDur - 1/Fe).';
+testSig = sin(2*pi*440*t);
+yDelay = effectDelay(testSig, tauSec, g, Fe);
+soundsc(yDelay, Fe); pause(testDur+0.5);
+
+%% ========================================================================
+%% Q 3.17 – Q 3.18 — Delay boucle + filtre moyenne K
+K = 10;
+yDelayFilt = effectDelayFiltre(testSig, tauSec, g, K, Fe);
+audiowrite('pianoDelayFiltre.wav', yDelayFilt, Fe);
+soundsc(yDelayFilt, Fe);
+disp('Q 3.18 Fichier pianoDelayFiltre.wav genere');
+
+%% ========================================================================
+%% Q 3.19 — Attenuation des hautes frequences (boucle filtree)
+w = linspace(0, pi, 2048);
+Hf = (1 - g * exp(-1j*w*tauEch)).^(-1) .* (sin(K*w/2)./(K*sin(w/2)));
+figure('Name','Q3-19 Attenuation HF');
+plot(w/pi*Fe/2, 20*log10(abs(Hf))); grid on
+xlabel('f (Hz)'); ylabel('dB');
+title('Boucle avec moyenne glissante — aigus attenues');
+
+%% ========================================================================
+%% === Fonctions locales (sans underscore) ================================
+function y = effectReverb(x,h)
+ y = filter(h,1,x); % convolution lineaire directe
+end
+
+function y = effectReverbFFT(x,h)
+ L = 2^nextpow2(length(x)+length(h)-1);
+ y = real(ifft( fft(x,L) .* fft(h,L) ));
+ y = y(1:length(x));
+end
+
+function h = impulseDelay(g,tauEch,Nmax)
+ n = 0:Nmax;
+ h = zeros(1,Nmax*tauEch+1);
+ h(n*tauEch+1) = (-g).^n;
+end
+
+function [b,a] = coeffDelay(g,tauEch)
+ b = [1 zeros(1,tauEch)]; % numerateur
+ a = [1 zeros(1,tauEch-1) -g]; % denominateur
+end
+
+function y = effectDelay(x,tauSec,g,Fe)
+ tauEch = round(tauSec*Fe);
+ [b,a] = coeffDelay(g,tauEch);
+ y = filter(b,a,x);
+end
+
+function y = effectDelayFiltre(x,tauSec,g,K,Fe)
+ tauEch = round(tauSec*Fe);
+ buf = zeros(tauEch,1);
+ y = zeros(size(x));
+ coef = ones(1,K)/K; % moyenne glissante
+ for n = 1:length(x)
+ delayed = buf(1);
+ buf = [buf(2:end); 0];
+ feedback = filter(coef,1,g*delayed);
+ y(n) = x(n) + feedback(1);
+ buf(end) = y(n);
+ end
+end
+
+function plotSpectrum(x,Fe,color,label)
+ N = length(x); f = (-N/2:N/2-1)*Fe/N;
+ X = fftshift(fft(x.*hann(N)));
+ SdB = 20*log10(abs(X)/max(abs(X))+eps);
+ plot(f,SdB,color,'LineWidth',1.2,'DisplayName',label);
+end
--
GitLab