Signalerx7


ØVELSE 7
i 31610


2009-04-13

Rasmus Berg Palm s062241


Mads Andersen s062237


Jacob Bjerring Olesen s062248


Indholdsfortegnelse

edit sections

Beskeder
til jacob i 2eren: Tidsforskellen bestemmes ikke som maksimumværdien af krydskorrelationen, men ved det korresponderende lag!
nej maks af krydskorrelationen er ikke givet i lag. og lag er tid
sigx7
Here is a place for your title Click me to edit ! false

1

load signals_white

fs = 16e3; % [hz] sampling frequency

%soundsc(recorded, fs)
x1 = recorded; %Data containing both speech and noise
x2 = noise; %Noise signal

2

(1)
\begin{aligned} x_{1}(n) &= s(n)*h(n) + g(n) \\ x_{2}(n) &= s(n) \\ \end{aligned}

Ved at beregne krydskorrelationen mellem to signaler fås et udtryk for middelsammenhængen mellem signalerne, hvor det ene i dette tilfælde er tidsforskudt i forhold til det andet. Størrelsen af tidsforskellen kan bestemmes ud fra laget der korresponderer med maksimumværdien af korrelationen. Korrelationen er de to signaler er findes ved følgende udtryk. Det benyttes at $r_{xy} &= x(l) * y(- l)$, samt at $r_{xx} &= x(l) * x(- l)$

(2)
\begin{aligned} r_{x1x2} &= x_{1}(l) * x_{2}(- l) \\ r_{x1x2} &= ( s(l)*h(l) + g(l)) * s(-l) \\ r_{x1x2} &= r_{ss}*h(l)+r_{gs}\\ \end{aligned}

Fra MATLAB-scriptet er maksimum fundet til at ligge efter 0.0059 sekunder, hvilket også er vist grafisk på figuren.

corr3.png

Dette resultat stemmer fint overens med den teoretiske tidsforskydning, der blev beregnet ved at dele afstanden mellem mikrofonerne (2 meter) med lydens hastighed (340 m/s), som ligeledes gav 0.0059 s. (se figur)

fs = 16e3; %hz

%theoretical delay
dist = 2; %m
v = 340.2; %m/s
tdelay = dist/v %s

%Calculating the cross-corelation of the two signals. 
[Rx1x2, k] = xcorr(x1, x2);% k is a vector of lag indices

%lags in seconds
tlags = k/fs;

%the calculated delay:
[m, delayi] = max(Rx1x2);% Finding the indices of the maximum value
delayt = tlags(delayi)% Calculated delay given in sec.

plot(tlags,Rx1x2)
xlim([4e-3 7e-3])
hold on
h = line([tdelay tdelay],[-1000 5000]); set(h,'Color','r');%Teoretisk delay
hold off
legend(['Correlations signal'],['Theoretically delay'])
xlabel('time [s]')
ylabel('r_{x1x2}')
title('Estimation of displacement')

3

For at få den bedste annullering af det målte støj signal, $x_{2}(n)$, skal det føres igennem et filter $w(n)$ givet ved det samme impuls respons som $h(n)$. Dette kan vises matematisk som følger…

(3)
\begin{aligned} x_{1}(n) - x_{2}(n) * w(n) &= y(n) \\ s(n)*h(n) + g(n) - s(n)*w(n) &= y(n)\\ \end{aligned}

Den bedste annullering af støjen ses når $y(n) &= g(n)$, da outputtet kun vil indeholde det talte signal, hvilket vil forekomme når $h(n) &= w(n)$.

4

Udtryk for $\Gamma_{yy}$ som funktion af W(f), H(f), $\Gamma_{gg},\Gamma_{ss}$:

Det udnyttes her de kommutative, assiciative og distributive egenskaber ved foldning, samt at $r_{xy}(l)= x(l)*y(-l)$.

(4)
\begin{aligned} y(n) &= s(n)*(h(n)-w(n))+ g(n) \quad \Rightarrow\\[3mm] r_{yy}(l) &= y(l)*y(-l) \\ &= [s(l)*(h(l)-w(l))+ g(l)]*[s(-l)*(h(-l)-w(-l))+ g(-l)]\\ &=r_{ss}(l)*[h(l)-w(l)]*[h(-l)-w(-l)]\\ &\quad +r_{sg}(l)*(h(l)-w(l))\\ &\quad +r_{gs}(l)*(h(-l)-w(-l))\\ &\quad +r_{gg}(l) \end{aligned}

Da s(n) er hvid, tilfældig støj er s(n) og g(n) totalt ubeslægtede, og krydskorrelationerne $r_{sg}(l), r_{gs}(l)$ kan derfor antages at have ubetydeligt små værdier. Dermed forsvinder de to midterste led. $\Gamma_{yy}(f)$ fås ved fouriertransformation af $r_{ss}(l)$:

(5)
\begin{aligned} r_{yy}(l) &=r_{ss}(l)*[h(l)-w(l)]*[h(-l)-w(-l)] + +r_{gg}(l) \quad \leftrightarrow\\[3mm] \Gamma_{yy}(f) &= \Gamma_{ss}(f)[H(f)-W(f)][H^*(f)-W^*(f)]+\Gamma_{gg}(f)\\ &= \Gamma_{ss}(f)[|H(f)|^2+|W(f)|^2-H^*(f)W(f)-W^*(f)H(f)]+\Gamma_{gg}(f) \end{aligned}

Her er bla. udnyttet at $h(-n) \leftrightarrow H(-f)=H^*(f)$ for et reelt signal h(n), og at foldning i tidsdomænet svarer til multiplikation i frekvensdomænet.

Indsættelse af h(n) og w(n):

(6)
\begin{aligned} h(n)&=\alpha \cdot \delta(n-d)\\ w(n)&=\beta \cdot \delta(n-k)\\[3mm] &\quad \leftrightarrow\\[3mm] H(f)&=\sum_{n=-\infty}^{\infty}\alpha \cdot \delta(n-d)e^{-j2\pi fn}=\alpha e^{-j2\pi fd}\\ W(f)&=\sum_{n=-\infty}^{\infty}\beta \cdot \delta(n-k)e^{-j2\pi fn}=\beta e^{-j2\pi fk}\\ \end{aligned}

Ved indsættelse i udtrykket for $\Gamma_{yy}(f)$ fås:

(7)
\begin{aligned} \Gamma_{yy}(f) &= \Gamma_{ss}(f)[\alpha^2+\beta^2-\alpha e^{j2\pi fd}\cdot\beta e^{-j2\pi fk}-\beta e^{j2\pi fk}\cdot \alpha e^{-j2\pi fd}]+\Gamma_{gg}(f)\\ &=\Gamma_{ss}(f)[\alpha^2+\beta^2-\alpha\beta (e^{j2\pi f(d-k)}+e^{-j2\pi f(d-k)})]+\Gamma_{gg}(f)\\ \end{aligned}

Ved brug af Eulers formel for cos fås:

(8)
\begin{align} \Gamma_{yy}(f) =\Gamma_{ss}(f)[\alpha^2+\beta^2-2\alpha\beta \cos(j2\pi f(d-k))]+\Gamma_{gg}(f) \end{align}

5

y(n) udregnes ud fra x1 og x2 som:

(9)
\begin{aligned} y(n) &=x_1(n)-x_2(n)*w(n)\\ &=x_1(n)-x_2(n)*(\beta \cdot \delta (n-k))\\ &=x_1(n)-\beta \cdot x_2(n-k) \end{aligned}

Ved at udføre den simple foldning $x_2(n)*w(n)$ ses at filtreringen svarer til at forsinke støjsignalet og skalere det. Dette passer fint med situationen der ønskes modelleret idet støjen ved mikrofonen højst sandsynligt er en dæmpet og forsinket udgave af støjen ved støjkilden. Dermed bør beta også være mindre end eller lig 1.
At forsinke og skalere et signal er mindre krævende for matlab end at folde to signaler, derfor ser find_y ud som følger:
function y =find_y(x1, x2, beta, k)

function y =find_y(x1, x2, beta, k)

%The conv0lution of w and x2 equals delaying x2 by k, and multiplying by
%beta:
x2convw = beta*[zeros(k-1,1); x2];

%To be able to substract x2convw from y, they must have equal lengths:
x2convw = x2convw(1:length(x1));

% y
y = x1 - x2convw;

Ved at bruge find_y med k=d og for forskellige beta fik vi følgende:

ynbeta01.png
ynbeta05.png
ynbeta075.png
ynbeta1.png

Her ses at med beta=0.75 fås et signal der ser ud til at have meget lidt støj. Vi lyttede til signalerne for at tjekke den påstand, og det var ganske rigtigt at beta=0.75 gav et meget fint lydsignal. Med beta=alpha er der perfekt kancellering af støjen. Vi estimerer dermed alpha=0.75

Følgende matlab kode blev brugt:

% ex 5
%lag in samples, in the signal. + 1 because t=0 correspondons to sample 1
d = delayt*fs+1

%Some different beta:
beta = [0.1 0.5 0.75 1];

%The corresponding y(n), we plot, and listen to:
for i=1:4
y = find_y(x1, x2, beta(i), d);
figure;
plot(y)
xlabel('sample number')
ylabel('y(n)')
title({['y(n) with beta=' num2str(beta(i))]})
soundsc(y,fs)
pause
end

6

Fra opg 2 er $r_{x_1x_2}=r_{ss}*h(n)+r_{gs}(n)$ hvilket i frekvensdomænet giver: $\Gamma_{x_1x_2}=\Gamma_{ss}H(f)+\Gamma_{gs}(f)$. Som omtalt i opg. 4 er værdierne $r_{gs}(l) \approx 0$, dermed bliver værdierne $\Gamma_{gs}(f) \approx 0$, da $\sum_{l=-\infty}^{\infty} 0 \cdot e^{-j2\pi fl}=0$. Da $x_2(n)=s(n)$ fås altså:

(10)
\begin{align} W_{opt}(f)=\frac{\Gamma_{x_1x_2}(f)}{\Gamma_{x_2x_2}(f)}=\frac{\Gamma_{ss}(f)H(f)}{\Gamma_{ss}(f)}=H(f) \end{align}

7

Følgende kode blev benyttet i denne opgave:

%% Ex 7
clear all; close all; clc;
load signals_white
fs = 16e3; %hz
x1 = recorded;
x2 = noise;
[Px1x2, W] = cpsd(x1,x2,[],512,1024, fs, 'twosided'); %Fourier transform of crosscorrelation of x1 and x2
[Px2x2, W] = cpsd(x2,x2,[],512,1024, fs, 'twosided'); %Fourier transform of autocorrelation of x2
Wopt = Px1x2./Px2x2; %Optimal filter spectrum
w = ifft(Wopt); %Optimal filter impulse response

k = conv(x2,w); %Reconstructed noise
y = x1 - k(1:length(x1)); %Reconstructed clean signal by subtracting reconstructed noise

wt = zeros(1,1024); %Ideal version of the filter w(t)
wt(95)=0.7518; %Values, 95 and 0.7518 read from a plot of w(t)
kt = conv(x2,wt); %Ideal reconstructed noise
yt = x1 - kt(1:length(x1)); %Ideal reconstructed signal

figure %Plotting
subplot(2,2,1)
plot(w)
xlabel('discrete time')
ylabel('Amplitude [V]')
title('Estimated w(t) by ifft(PSD(x1x2)/PSD(x2x2))')
ylim([-0.1 1])
subplot(2,2,2)
plot(y);
xlabel('discrete time')
ylabel('Amplitude [V]')
title('x1(t) filtered with w(t)')

ylim([-1 1])
subplot(2,2,3)
plot(wt)
xlabel('discrete time')
ylabel('Amplitude [V]')
title('wt(t), Ideal form of w(t): beta = 0.7518, k = 0.0059s')
ylim([-0.1 1])
subplot(2,2,4)
plot(yt);
xlabel('discrete time')
ylabel('Amplitude [V]')
ylim([-1 1])
title('x1(t) filtered with wt(t)')

SNRpre = (sqrt(mean(speech.^2))/sqrt(mean((x1-speech).^2)))^2 %SNR pre-filtering
SNRpost = (sqrt(mean(speech.^2))/sqrt(mean((y-speech).^2)))^2 %SNR after filtering with w(t)
SNRpost_t = (sqrt(mean(speech.^2))/sqrt(mean((yt-speech).^2)))^2 %SNR after filtering with ideal w(t)

Hvilket gav følgende output:

ex7.png
SNRpre =

    0.4524

SNRpost =

   10.8145

SNRpost_t =

  344.5720

Forsinkelsen, k, blev fundet som den diskrete tid for hvilke w(t)'s impulsrespons havde sin maksimale værdi og alpha blev fundet som denne maksimale værdi. Den ideelle form af w(t) er jo netop en impulsfunktion med k forsinkelse og alpha amplitude.
Signal støj forholdet (SNR) taler sit tydelige sprog, før filtrering var der et SNR på 0.45, altså over halvdelen af signalets energi var støj, hvilket gjorde det umuligt at høre hvad der blev sagt. Efter filtrering var SNR på 10.8, hvilket gjorde at man tydeligt kunne høre hvad der blev sagt. Ved filtering med det ideelle filter var signal støj forholdet imponerende 344.6 hvilket ikke kunne skelnes fra det rene signal (speech vektoren) med det blotte øre.

8

Følgende kode blev benyttet til denne opgave:

%% Ex 8
%Mostly a copy of Ex 7. See Ex 7 for additional comments.
clear all; close all; clc;
load signals_hairdryer
% fs = 16e3; %hz
x1 = recorded;
x2 = noise;
[Px1x2, W] = cpsd(x1,x2,[],512,1024, fs, 'twosided');
[Px2x2, W] = cpsd(x2,x2,[],512,1024, fs, 'twosided');

Wopt = (Px1x2)./Px2x2;
w = ifft(Wopt);

k = conv(x2,w);
y = x1 - k(1:length(x1));
wt = zeros(1,1024);
wt(95)=0.7518;
kt = conv(x2,wt);
yt = x1 - kt(1:length(x1));

figure %Plotting
subplot(2,2,1)
plot(w)
xlabel('discrete time')
ylabel('Amplitude [V]')
title('Estimated w(t) by ifft(PSD(x1, x2)/PSD(x2, x2)))
ylim([-0.5 1])
subplot(2,2,2)
plot(y);
xlabel('discrete time')
ylabel('Amplitude [V]')
title('x1(t) filtered with w(t)')

ylim([-1 1])
subplot(2,2,3)
plot(wt)
xlabel('discrete time')
ylabel('Amplitude [V]')
title('wt(t), Ideal form of w(t): beta = 0.7518, k = 0.0059s')
ylim([-0.5 1])
subplot(2,2,4)
plot(yt);
xlabel('discrete time')
ylabel('Amplitude [V]')
ylim([-1 1])
title('x1(t) filtered with wt(t)')

SNRpre = (sqrt(mean(speech.^2))/sqrt(mean((x1-speech).^2)))^2
SNRpost = (sqrt(mean(speech.^2))/sqrt(mean((y-speech).^2)))^2
SNRpost_t = (sqrt(mean(speech.^2))/sqrt(mean((yt-speech).^2)))^2

Hvilket gav følgende output:

ex8.png
SNRpre =

   12.7585

SNRpost =

    5.7308

SNRpost_t =

  8.3061e+003

Det ses altså at filtreringen med filtered fundet ved w=ifft(PSD(x1, x2)/PSD(x2, x2)) ikke hjælper. Tværtimod ser det ud til at SNR falder. Hvis man lytter til signalerne er der ikke den store forskel dog. Ved filtrering med det ideelle filter fra excercise 7 fås dog et meget imponerende SNR på 8306. Man skulle tro at det bedste ideelle filter ville findes ved aflæsning på den nye w(t), men det ses på figuren at denne ikke ligner et impulsrespons særligt meget. Der må være noget galt. Ved gennemgang af teorien ses det at der er en central forudsætning for at Wopt(f)=PSD(x1, x2)/PSD(x2, x2), netop at støjsignalet og talesignalet er fulstændigt ukorreleret. Ved tilføjelse af disse linjer kode:
[Pgs, W] = cpsd(speech,x2,[],512,1024, fs, 'twosided');

Wopt = (Px1x2-Pgs)./Px2x2;

Trækkes PSD af støjsignalet og talesignalet fra PSD af x1 og x2 inden der divideres med PSD(x2,x2). Altså burde man få et mere præcist udtryk fra H(f) hvis at PSD(speech,x2) ikke kan antages at være ubetydeligt lille. Dette giver:
SNRpre =

   12.7585

SNRpost =

  6.9350e+003

SNRpost_t =

  8.3061e+003

Altså en imponerende forbedring i virkningen af filteret w(t). Det må konkluderes at den grundlæggende forskel var at talesignalet og støjsignalet på den ene eller den anden måde var korreleret.
I en realistisk målesituation er man desværre ikke så heldig at have støjsignalet, signalet og det kombinerede signal, hvorfor denne metode ikke kan benyttes ved enhver lejlighed til at opnå fantastiske signal-støj forhold.

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License