Upload New File

PARTIE3.m

+%% ========================================================================
+%  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