Signalerx3


ØVELSE 3
i 31610


2009-03-16

Rasmus Berg Palm s062241


Mads Andersen s062237


Jacob Bjerring Olesen s062248


Indholdsfortegnelse

Teori forberedelse

1

En ligefordeling har sandsynlighedstæthedsfunktionen:

(1)
\begin{align} w_X(\xi)=\frac{1}{b-a}, \quad for \quad a<\xi <b, \quad \quad 0, \quad ellers \end{align}

Hvor X er en stokastisk variabel taget fra ligefordelingen. Den forventede værdi af X, eller middelværdien for ligefordelingen, kan findes ved formel (8.13) i noterne:

(2)
\begin{align} E\{ X\} =\int_{-\infty}^\infty \xi w_X(\xi ) d\xi \end{align}
(3)
\begin{align} =\int_{a}^b \frac{\xi }{b-a}d\xi \end{align}
(4)
\begin{align} =\frac{b+a}{2} \end{align}

Variansen fås af (8.16):

(5)
\begin{align} \sigma^2\{ X\} =\int_{-\infty}^\infty (\xi - E\{ X\} )w_X(\xi ) d\xi \end{align}
(6)
\begin{align} =\int_{a}^b (\xi - \frac{b+a}{2})\frac{1}{b-a} d\xi \end{align}
(7)
\begin{align} =\frac{(b-a)^2}{12} \end{align}

Spredningen er derfor:

(8)
\begin{align} \sigma\{ X\} = \frac{b-a}{\sqrt{12}} \end{align}

2

Den forventede værdi for en stokastisk process $X(n)$, eller middelværdien af den fordeling, $w_X(\xi ; n)$ som $X(n)$ er samplet fra er:

(9)
\begin{align} E\{ X(n)\} =\int^\infty_{-\infty}\xi w_X(\xi; n)d\xi \end{align}

Autokorrelationen af en stationær, stokastisk process er:

(10)
\begin{align} R_X(k)=\int ^\infty _{-\infty}\int ^\infty _{-\infty}\xi _1 \xi _2 w_X(\xi _1, \xi _2; n, n+k)d\xi _1 d\xi _2 \end{align}

Det ses at autokorrelationen af et stokastisk signal fra en stationær, stokastisk process er middelværdien af produktprocessen $X(n)X(n+k)$ fra den todimensionelle sandsynlighedsfordeling med sandsynlighedstæthedsfunktionen $w_X(\xi _1, \xi _2, n, n+k)$, dvs:

(11)
\begin{align} R_X(k)=E\{X(n)X(n+k)\} \end{align}

heraf følger det at

(12)
\begin{align} R_X(0)=E\{X^2(n)\} \end{align}

Variansen af en stokastisk process er

(13)
\begin{align} \sigma^2\{ X(n) \}=E\{ (X(n)-E\{ X(n)\} )^2\} = E\{X^2(n)\} - E^2\{ X(t)\} \end{align}

dvs

(14)
\begin{align} \sigma^2\{ X(n) \} = R_X(0)- E^2\{ X(t)\} \end{align}

eller

(15)
\begin{align} R_X(0)=\sigma^2\{ X(n) \} + E^2\{ X(t)\} \end{align}

Sandsynlighedstæthedsfunktionen for en normalfordeling med middelværdi m og spredning s, er givet ved:

(16)
\begin{align} w_X(\xi ) = \frac{1}{\sqrt{2\pi s^2}}\exp{(-\frac{(\xi - m)^2}{2s^2})} \end{align}

autokorrelationen i k= 0 af en process derfra er altså specifikt:

(17)
\begin{equation} R_X(0)=s^2 + m^2 \end{equation}

3

Den teoretiske autocorrelation for

(18)
\begin{equation} Y(n) = r_{h}(n) * (X(n)-E(X(n)) \end{equation}

skal findes.
hvor $r_{h}(n)$ er givet ved:

(19)
\begin{align} r_{h}(n) = rect(n/20) = \frac{1}{20}, \quad for \quad 0 < n \leq 20, \quad \quad 0, \quad ellers \end{align}

For fremtiden vil $X(n)$ betegne $X(n)-E(X(n))$
Ifølge 8.168 gælder da følgende:

(20)
\begin{equation} R_{Y}(k) = R_{h}(k) * R_{X}(k) \end{equation}

8.179 giver autocorrelationfunktionen for digital hvid støj, $X(n)$

(21)
\begin{align} R_{X}(k) = A\delta(k) \end{align}

A kan her findes, da:

(22)
\begin{align} A = P_{X} = R_{X}(0) = E(X(n))^2 + \sigma(X(n))^2 \end{align}

Og ved indsættelse af værdier fås:

(23)
\begin{align} E(X(n)) = 0, \sigma(X(n))^2 = \frac{(b-a)^2}{12}, b = 1/2, a = -1/2 \end{align}
(24)
\begin{align} R_{X}(k) = \frac{1}{12}\delta(k) \end{align}

Hvor at a og b er grænserne for X's PDF. De var 0 og 1, men ved at trække middelværdien (1/2) fra blev de -1/2 og 1/2. Middelværdien er indlysende lig 0.
$R_{h}(k)$ findes vha. 8.208:

(25)
\begin{align} R_{h}(k) = \frac{1-\frac{|k|}{20}}{20}, \quad for \quad |k|<20, \quad \quad 0, \quad ellers \end{align}

Nu er blot tilbage at folde $R_{h}(k)$ og $R_{X}(k)$. Da $R_{X}(k)$ er en delta funktion med amplituden 1/12, giver dette blot:

(26)
\begin{align} R_{Y}(k) = \frac{1}{12}\frac{1-\frac{|k|}{20}}{20}, \quad for \quad |k|<20, \quad \quad 0, \quad ellers \end{align}

Opgaver

1

For at estimere en sandsynlighedstæthedsfordelingen for et diskret signals amplitude, inddeles amplitudevidden (altså det område som amplituden i signalet spænder over) i et passende antal intervaller, og det tælles hvor mange signalværdier der falder inde for hvert interval. Denne optælling kan gøres i matlab vha. funktionen histc. Den relative hyppighed af et tal fra hvert interval fås ved at dividere de optællede hyppigheder med det samlede antal af værdier. For at få tætheden divideres vidden af hvert interval. I funktionen probdens.m er dette gjort, hvor antallet af intervaller er 100:

function [pd, iv] = probdens(sig, doplot)

% The greatest value in the signal:
maks = max(sig)

% The smallest value in the signal:
mini = min(sig)

% The edges of the intervals:
edges = linspace(mini, maks, 101); 

counts = histc(sig, edges); % Counts the values in the interval [edges(i);edges(i+1)[, 
%for i=1:100, and place the counts in af 100 colum vector

% The next is because the the last element in counts using histc is the 
% number of values in sig that matches edges(end) - we want these values 
% counted in the interval [edges(end-1);edges(end)] instead:
counts(end-1) = counts(end-1)+counts(end); 
counts = counts(1:end-1);

% The occurence of the different values:
occur = counts/length(sig);

% The width of each interval:
width = (maks-mini)/(length(edges)-1);

% The probability density::
pd = occur/width;

% amplitude axis (interval values):
iv = width/2:width:maks-width/2;

if nargin > 1
    figure;
    bar(iv, pd)
    xlabel('signal amplitude')
    ylabel('occurence')
end

Med funktionen rand fra MATLAB er lavet et tilfældigt signal X med 10000 værdier. Signalet er plottet i forhold til tiden, idet F_s=1000Hz. Plottet ses nedenfor, og det kan ses at signalet består af en masse værdier mellem 0 og 1. Dette skyldes at rand netop udtager værdier fra en ligefordeling mellem 0 og 1.

opg1xt.png

En ligefordeling mellem 0 og 1 har sandsynlighedstæthedsfunktionen $1 for 0<\xi <1, \quad 0, \quad ellers$, jf. teoriafsnit 1.
Sandsynlighedstætheden for signalets amplitude er estimeres vha. probdens.m og nedenfor ses denne plottet sammen med sandsynlighedstæthedsfunktionen for en ligefordeling.

opg1ligepd.png

De er ikke helt ens da der er tilfældighed involveret og der trækkes et endeligt antal samples - i princippet kunne rand() trække det samme tal igen og igen, hvilket ville øge forskellen, men sandsynligheden for at det sker er lille. Generelt set vil de nærme sig hinanden jo større mængder samples der udtages.
Nedenfor ses dette plottet:

ex1errprsamples.png

Følgende matlab kode blev brugt til at producere dette plot:

clear all; close all; clc;
figure
hold on
for i=10:10:1000
  sig = rand(i,1);
  aerr = mean(abs(1-probdens(sig)));
  plot(i,aerr,'.')
end
hold off
xlabel('Samples'); 
ylabel('Mean abs diff'); 
title('Mean abs diff between analytical pdf and computed pdf as a function of sample size');

2

Et lavpass, lineært filter med rektangulært impulsrespons og forstærkning på 1 ved dc-komponenten, fås med:

(27)
\begin{align} h(n) = 1/N, \quad \leq n \leq N-1 \leftrightarrow \end{align}
(28)
\begin{align} H(F)=\frac{1}{N}\frac{\sin(\pi F/F_sN)}{\sin(\pi F/F_s)}e^{-j\pi F/F_s(N-1)} \end{align}

At impulsresponset er rektangulært ses direkte. At forstærkningen er 1 ved dc-komponenten kan ikke ses ved direkte indsættelse af F=0 i H(F), men ved at taylorudvikling af 1. orden af tæller og nævner :

(29)
\begin{align} H(0) = \lim_{F \rightarrow 0}\frac{1}{N} \frac{\sin (0)+\pi /F_sN \cos (0)F}{\sin (0)+\pi /F_s \cos (0)F}=\lim_{F \rightarrow 0}\frac{1}{N}N=1 \end{align}

H(F) er 0 første gang når:

(30)
\begin{align} \pi F/F_sN=\pi \quad \Leftrightarrow \quad N=F_s/F \end{align}

Så hvis det første 0 skal være ved $F=50$Hz og $F_s=1000$Hz, skal

(31)
\begin{equation} N=1000Hz/50Hz=20 \end{equation}

Filterets overføringsfunktion ses plottet nedenfor.

opg2HF.png

Signalet $X(n)$ fra opgave 1 er fra en ligefordeling med $0 <\xi <1$. Middelværdien for denne fordeling er (jf. teoriafsnit 1)

(32)
\begin{align} E\{X(n)\}=1/2 \end{align}

Signalet $X_2(n) = X(n) -1/2$ er blevet filtreret med det rektangulære filter, ved en foldning af firkantfunktionen og signalet. Det filtrerede signal $Y(n)$ ses sammen med $X(n)$ nedenfor.

opg2XY.png
opg2XYzoom.png

Det ses at Y flukturerer omkring 0, hvilket er klart nok eftersom middelværdien blev trukket fra, og dc-komponenten af signalet passere uhindret gennem filteret. Derudover ses det at amplituden af Y generelt er mindre end amplituden af X. Dette passer fint med at kun dc-komponent passere uhindret gennem filtret og resten af frekvenserne dæmpes meget kraftigt, jf. figuren af H(F). På zoom figuren ses at de høje frekvenser (hurtige fluktationer) er meget mere dæmpede end de lave frekvenser (langsomme fluktationer), hvilket også passer med filterets overføringskarakteristik.

Fordelingen af det filtrerede signal Y er anskueliggjort vha. probdens.m fra opg. 1 og ses på figuren nedenfor.

opg2pd.png

Det ses at elementerne i Y tilnærmelsesvis er normalfordelt. Dette skyldes at der foldes med en firkantfunktion. X2 foldet med en firkantfunktion med længden 20, giver en funktion Y, hvis elementer hver er en sum af 20 elementer i X2. X2 består af uafhængige stikprøver fra en uniform fordeling middelværdi $\mu=0$ og varians $\sigma^2=4/12=1/3$ (jf. teori afsnit 1). Hvert element i Y er således en sum af 20 samples fra en uniform fordeling med $\mu=0$ og varians $\sigma^2=4/12=1/3$. Fra den centrale grænseværdi sætning (fra statistik og sandsynlighedsregning) haves at udtages stikprøver fra en hvilken som helst fordeling, som hver især summeres, vil summene fordele sig efter en normalfordeling, når stikprøvestørrelsen går mod uendelig. Her er stikprøvestørrelsen 20, og derfor ses elementerne i Y ikke perfekt normalfordelt, men approksimationen ses.
Normalfordelingen der konvergeres mod, har middelværdi som fordelingen der samples fra (dvs 0) og spredning $\sigma /\sqrt(20)$, og denne fordeling er også indtegnet på figuren (se s 212 i Miller and Freunds Probability afnd Statistics for Engineers).

Følgende MATLAB-kode er brugt til opg 1 og 2:

% Ex1
clc; clear all; close all;

% 10 000 samples from a uniform probability distribution:
X = rand(1e4,1); 

% Sampling frequency:
fs = 1000; %Hz

% Time axis:
t = 0:1/fs:(length(X)-1)/fs; 

% Plot of signal:
figure;
plot(t,X)
xlabel('Time (s)')
ylabel('amplitude')
title('The random signal X(t)')

% The probability density of the signal amplitude:
[dens,iv] = probdens(X);

% Theoretical probability density function
thedens = 1; 

figure;
bar(iv,dens)
hold on
plot([0 1],[thedens thedens], 'k-','LineWidth', 2)
xlabel('signal amplitude')
ylabel('occurence')
legend('Pd of signal', 'Theoretical pd')
title('Probability densities')
hold off

% Ex2

% The impulse response of the lowpass filter:
N = 20;
h(1:N)= 1/N;

% Plot of the amplitude spectrum of the transfer function:
addpath('C:\Users\Mads\Desktop\mads\mads\dokumenter\med og tek\anvendt signalbehandling\ex2\matlab')
[H, F] = plot_spectrum(h,fs, 1000);
figure;
plot(F,H)
legend('H(F)')
xlabel('F (Hz))')
ylabel('amplitude')
title('The transfer function of the filter')

% The filtering:
Y = conv(h,X - 1/2);

% Timeaxis for the filtered signal:
tY = 0:1/fs:(length(Y)-1)/fs;

% Plot of unfiltered and filtered signals
figure;
plot(t,X, 'b-')
hold on
plot(tY,Y, 'r-')
xlabel('Time (s)')
ylabel('amplitude')
legend('X','Y')
title('The raw signal (X) and the filtered signal (Y)')
axis([0 10 -0.3 1.2])
hold off

% pd:
[dens,iv] = probdens(Y);

% Theoritical pdf:
mu    = 0;
sigma = (4/12)/sqrt(20);
gauss = normpdf(iv, mu, sigma);

% Figure of the to pd's:
figure;
bar(iv, dens)
hold on
plot(iv, gauss, 'r', 'LineWidth', 2)
xlabel('signal amplitude')
ylabel('occurence')
legend('Pd of signal', 'Theoretical pdf')
title('Probability densities')
hold off

3

ex3corr.png

Ovenfor ses en sammenligning af den teoretiske og den udregnede autocorrelation af Y(n), som defineret i teoriafsnit 3.
Der ses god sammenlignelighed, dog med visse afvigelser. Afvigelserne er et resultat af at den teoretiske udregnede correlation regner på sandsynelighedsfordelingen for den stokastiske process, hvorimod den komputationelle regner på et endeligt antal samples fra fordelingen. Den gennemsnitlige absolutte forskel mellem disse 2 er proportionel med antal samples udtaget, som vist før.
Følgende matlab kode blev brugt til ovenstående:

clear all; close all; clc;
ex1
close all; clc;
[Ryc, lags] = xcorr(y,100,'unbiased');
i = 1;
for k=lags
    if(abs(k)<20)
        Ryt(i) = 1/12.*(1-(abs(k)./20))./20;
    else
        Ryt(i)  = 0;
    end
    i = i+1;
end
plot(lags, Ryt, lags, Ryc);
xlabel('Lag (k)')
ylabel('Correlation between Ry(n) and Ry(n+k)');
title('Comparison of theoretical and computational correlation of Ry(n) and Ry(n+k)');
legend('Theoretical','Computational');

4

Ultralydssignalets sandsynlighedstæthedsfunktion er sammenlignet med en normalfordeling med middelværdi og varians udregnet fra ultralydssignalet. Resultatet ses nedenfor.

opg4.png

Det ses at ultralydssignalet tilnærmelsesvist er normalfordelt. Dette har at gøre med at hvert element i signalet stammer fra utrolig mange små bidrag der summeres. Selvom at disse bidrag ikke er fra samme fordeling, vil summen af dem stadig (under visse antagelser) konvergere mod en normalfordeling ifølge den centrale grænseværdisætning.

Følgende matlabkode er brugt til opg 4:

% Ex4
close all; clear all; clc;
load ult_sig2

% The signal plotted
figure;
plot(ult_sig2)

% The mean value of the signal:
mu  = mean(ult_sig2);
sigma = std(ult_sig2);

% The pd of the signal:
[dens,edges] = probdens(ult_sig2);

% The pdf of a normal distribution with mean = mu, std = sigma:
gauss = normpdf(edges, mu, sigma);

% Figure of the to pdf's:
figure;
bar(edges, dens)
hold on
plot(edges, gauss, 'r', 'LineWidth', 3)
xlabel('signal amplitude')
ylabel('occurence')
legend('Pd of signal', 'Pdf of gauss pdf with mean and std from signal')
title('Probability densities')
hold off

5

ex5single.png

Ovenfor ses g, som er en 6Hz sinus med amplituden 1 og fasen 0, kontamineret med gaussisk støj (2*randn). Sinuskomponenten er meget svær at se.

Nu udregnes 100 realisationer af g og deres middelværdi findes:

ex5mean.png

Sinus komponenten ses nu tydeligt, og støjen lader til at være forsvundet.
Dette gælder da

  • Støjen ikke har nogen korrelation med sinus komponenten
  • En realisation af støjen ikke har nogen korrelation med en anden realisation.

Altså vil middelværdien i alle punkter for støjen gå mod 0, mens at middelværdien af sinuskomponenten vil være værdien af sinuskomponenten.
Den samlede middelværdi i alle punkter vil altså gå mod sinuskomponentens værdi for punktet for mængden af realisationer gående mod uendelig.

Følgende matlab kode blev benyttet til ovenstående:

clear all; close all; clc;
fs=1000; f0 = 6;
g = sin(2*pi*f0*(0:2000)/fs) + 2*randn(1,2001);
figure
plot(g);
title('6Hz sine contaminated with gaussian noise');
xlabel('Discrete time, n');
ylabel('Amplitude');
for i=1:100
  g(i,:) = sin(2*pi*f0*(0:2000)/fs) + 2*randn(1,2001);
end
figure
plot(mean(g))
title('Mean of 100 6Hz sines contaminated with gaussian noise');
xlabel('Discrete time, n');
ylabel('Amplitude');

edit sections

Here is a place for your title
Click me to edit !
Drag me !
sigx3
Here is a place for your title Click me to edit ! false
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License