Signalerx2


ØVELSE 2
i 31610


2009-03-10

Rasmus Berg Palm s062241


Mads Andersen s062237


Jacob Bjerring Olesen s06248


Indholdsfortegnelse

Teori forberedelse

1

Har indført $F_0$ i stedet for $f_0$ for frekvensen af tidssignalet. Dermed kan $f$ bruges som den relative frekvens af det digitale signal, altså $f_0=F_0/F_s$, og generelt er $f=F/F_s$.

Et periodisk, diskret-tid signal med perioden N, kan repræsenteres som fourierrækken (jf. (4.2.7),(4.2.8) i Proakis bog):

(1)
\begin{align} x_p(n)=\sum_{k=0}^{N-1}c_ke^{j2\pi kn/N}, \-\infty <n< \infty \end{align}

hvor

(2)
\begin{align} c_k=1/N\sum_{n=0}^{N-1}x_p(n)e^{-j2\pi kn/N}, k=0,1,...,N-1 \label{ck} \end{align}

Fourierkoefficienterne $c_k$ beskriver signalets i frekvensdomænet hvor basis er komplekse eksponentialfunktioner, altså signalets komplekse spektrum. Hvis $c_k$ afbildes som funktion af $f$ findes $c_1$ ved $f=1/N$, $c_2$ ved $f=2/N$, $c_0$ ved $f=0$ osv. Afbildes $c_k$ som funktion af $F$ findes $c_1$ ved $F=F_s/N$, $c_2$ ved $f=2F_s/N$, $c_0$ ved $F=0$ osv. generelt $c_k$ ved $F=kF_s/N$. Det komplekse spectrum / Fourier koefficienterne er periodiske med N, dvs $c_k+pN=c_k$ hvor p er et heltal. Ses $c_k$ som funktion af $F$ er det en periode på $F=NF_s/N=F_s$.
Det følger af ovenstående at spektret af én kompleks eksponentialfunktion med amplituden 1, er én fourierkoefficient - med amplituden 1, som gentager sig med perioden $F_s$ hvis spektret afbildes som funktion af $F$. Specielt er en eksponentialfunktion på formen

(3)
\begin{align} g(n)=exp(j2\pi \frac{F_0}{F_s}n) \end{align}

en periodisk funktion der kan skrives som $x_p$-udtrykket, med $c_k=1$. Det ses så at $F_0/F_s=k/N \Leftrightarrow kF_s/N=F_0$. Fourierkoefficienten, $c_k$, til eksponentialfunktionen findes som funktion af $F$ ved $F=kF_s/N$, dvs netop ved $F=F0$. Da spektret er periodisk med $F_s$ er kan spektret for $g(n)$ udtrykkes ved den diskrete enheds-samplefunktion som $\delta(F-F0+pF_s), p \in Z$.

3

1

(4)
\begin{align} h(n)=1/N, \quad 0\leq n\leq N-1, \qquad h(n)=0, \quad ellers \leftrightarrow \end{align}
(5)
\begin{align} H(F)=\sum_{n=0}^{N-1}\frac{1}{N}e^{-j2\pi F/F_sn}=\frac{1}{N}\sum_{n=0}^{N-1}(e^{-j2\pi F/F_s})^n \end{align}

En kvotientrække, dvs:

(6)
\begin{align} H(F)=\frac{1}{N}\frac{1-e^{-j2\pi F/F_sN}}{1-e^{-j2\pi F/F_s}} \end{align}
(7)
\begin{align} =\frac{1}{N}\frac{(e^{j\pi F/F_sN}-e^{-j\pi F/F_sN})e^{-j\pi F/F_sN}}{(e^{-j\pi F/F_s}-e^{-j\pi F/F_s})e^{-j\pi F/F_s}} \end{align}
(8)
\begin{align} =\frac{1}{N}\frac{\sin(\pi F/F_sN)}{\sin(\pi F/F_s)}e^{-j\pi F/F_s(N-1)} \end{align}

2

(9)
\begin{align} h(n)=\frac{1}{2}\delta (n) - \frac{1}{2}\delta (n-10) \Leftrightarrow \end{align}
(10)
\begin{align} h(n+5)=\frac{1}{2}\delta(n+5)-\frac{1}{2}\delta(n-5) \end{align}

Udnyttes linearitet og reglen om forskydning af tidsaksens nulpunkt, fås:

(11)
\begin{align} H(F)=(\frac{1}{2}e^{-j2\pi F/F_s(-5)}-\frac{1}{2}e^{j2\pi F/F_s5})e^{j2\pi F/F_s5}=j\sin (10\pi F/F_s)e^{j10\pi F/F_s} \end{align}

Opgave 1

Periodiske signaler og deres spektre i MATLAB

Periodiske signaler er ikke mulige at repræsentere i MATLAB da de har uendelig længde. Men det er faktisk muligt vha. MATLAB at finde det korrekte komplekse spektrum for et periodisk signal, hvis en periode af signalet er repræsenteret.

DFT (fft) af én periode af et periodisk signal med perioden N, er ((7.1.20 i Proakis bog):

(12)
\begin{align} X(k)=\sum_{n=0}^{N-1}x(n)e^{-j2\pi kn/N}, k=0,1,...,N-1 \end{align}

hvor $x(n)$ er en periode af et digitalt periodisk signal $x_p(n)$.
Spektret af et periodisk digitalt signal $x_p(n)$ er givet ved fourierkoefficienterne $c_k$ som udtrykt i Teori forberedelse afsnit 1. Ses på det udtrykt bides mærke i at der kun summeres over én periode, dvs. $x_p(n)$ kan erstattes af $x(n)$, hvorved det ses at

(13)
\begin{align} c_k=1/N\cdot X(k) \end{align}

Et periodisk signals komplekse spektrum kan altså findes ved DFT af én periode af signalet og en skalering med $1/N$. Det kan yderligere vises (jeg har dog ikke en reference til et sted hvor det er gjort) at hvis man har et helt antal perioder af et periodisk signal, kan det komplekse spectrum findes eksakt ved DFT-spektret X_p(k) som:

(14)
\begin{align} c_k=1/(pN)\cdot X(k), p \in \{1,2,3,...\} \end{align}

hvor p er antallet af perioder.
En implementering af hvordan man finder spektret af et periodisk signal, repræsenteret ved et helt antal perioder, er gjort i plot_spectrum2.m:

function [ckamp,f] = plot_spectrum2(sigpper, fs, doplot)
% Calculates the absolute value of the complex spectrum of a periodic, 
% digital signal. If there is 3 input arguments (doplot has a value) a plot
% of the spectrum is made. sigpper is the periodic signal represented by a
% an integer number of periods.

% The DFT of an integer number of periods of the periodic signal:
X = fftshift(fft(sigpper)); 
% By only giving fft one input argument, it makes as many equidistant 
% samples of the fourier spectrum as there is elements in signal.
% The fftshift frequency-shifts the spectrum, so the zero-frequency 
% component is the center of the spectrum, and the spectrum therefore can 
% be visualised with the zero-frequency component in the middle of the 
% spectrum.  

% The conversion factor: 
pN = length(sigpper);

% The Fourier coefficients (One period of the complex spectrum):
ck = 1/pN*X;

% The amplitude spectrum (one period):
ckamp = abs(ck);

% The frequency values at which ck is located:
f = [0:pN-1]/pN*fs-fs/2;
% f is a vector which contains N elements equally spaced in the half-open
% interval, [-fs/2 ; fs/2[. 

if(nargin>2 && doplot==1)
  figure
  stem(f,ckamp)
  legend('spec')
  xlabel('Frequency (F) [Hz]'),ylabel('Amplitude [Voltage]')
end

Plot af g(n) med varighed 1 sekund

Ønsker at repræsentere det periodiske signal

(15)
\begin{align} g(n)=exp(j2\pi F_0n/F_s) \end{align}

i MATLAB. Det skal repræsenteres med det antal samples der svarer til 1 sekund, og $F_s=1000$Hz - altså 1000 samples. Spektret af $g(n)$ ønskes plottet for $F_0=230$Hz, $F_0=700$Hz og $F_0=1750$Hz, vha. plot_spectrum2.m. For at plot_spectrum2.m giver de korrekte spektre skal $g(n)$ være repræsenteret ved et helt antal perioder.
Fra s. 15 er den fundamentale periode af et digitalt signal det mindste heltal N, som opfylder $f_0N=k, k\in \{1,2,3...\}$. $k$ og $N$ er altså de mindste heltal der opfylder $k/N=f_0=F_0/F_s$ - eller sagt på en anden måde er $k/N$-brøken er den mest forkortede udgave af $F_0/F_s$-brøken, dvs $F_0=c1\cdot k, \quad F_s=c1\cdot N$, hvor $c1 \in \{1,2,3,..\}$ (og må have enheden Hz for at enhederne stemmer overens). Så hvis man repræsenterer det periodiske digitale signal med Fs-samples, er der repræsenteret et helt antal perioder af det digitale signal.
Her er der repræsenteret 1 sekund af signalerne, hvilket præcist svarer til $F_s$ samples, idet $F_s$ netop er antallet af samples pr sekund, dvs. alle de periodiske signaler ($g(n)$ er repræsenteret i MATLAB ved et helt antal perioder, og de korrekte spektre kan derfor udregnes. Dette er gjort vha. plot_spectrum2.m og resultaterne ses nedenfor.

opg1_230hz.png
opg1_700hz.png
opg1_1750hz.png

Som tidligere omtalt, er spektret af et digitalt signal periodisk med perioden $F_s$, og al information omkring et spektrum fås derved ved at se spektret i fx intervallet $[-F_s/2;F_s/2]$, hvilket er det interval der er brugt på figurerne.
Som omtalt i teorien er spektret af g(n) lig med $\delta(F-F0+pF_s), p \in Z$. Ser man derfor i et interval med længden $F_s$, som her, forventes altså kun én $\delta$ på hver figur - og det er netop hvad der ses. På spektret for $g(n)$ med $F_0=230$Hz ses funktionen $\delta(F-230)$ -altså en diskret $\delta$ i $F=F_0$. På spektret for $g(n)$ med $F_0=700$Hz ses en $\delta$ i $F=-300=F_0-1\cdot F_s$, så det er også som forventet idet $F_0>F_s$ og således ikke kan være med på figuren. Ligeledes ses på spektret for $g(n)$ med $F_0=1750$Hz en $\delta$ i $F=-250=F_0-1\cdot F_s$.

MATLAB-koden brugt til denne opgave ses nedenfor:

% excercise 1
close all; clear all; clc;
% The duration of g(n):
dur = 1; % s
% The sampling frequency:
Fs = 1000;% Hz
% The n-vector:
n = 0:Fs*dur-1; 

% Generation of g(n) for different values of F0, and plotting the spectrum 
% using plot_spectrum2:
F0 = [230 700 1750]; % Hz

for i=1:3
g1 = exp(j*2*pi*F0(i)*n/Fs);
[spec, f]=plot_spectrum2(g1, Fs);
figure;
stem(f,spec)
legend('spectrum of g')
xlabel('Frequency (F) [Hz]'),ylabel('Amplitude [Voltage]')
title(['The spectrum of g(n) with F0=' num2str(F0(i)) 'Hz'])
end

Opgave 2

Signalet g(n) er en sum af rene toner med frekvenser 100Hz, 150Hz og 300Hz samt en DC komponent.

(16)
\begin{align} g(n) = sin(2\pi f_{0} n/fs) + 0.2 sin(2\pi f_{1} n/fs) - 0.5 sin(2\pi f_{2} n/fs) + 1.2 \end{align}

Hvis den største fælles frekvens, Fb, findes vil den mindste fælles periode, Tb, være fundet da Tb = 1/Fb. Det ses umiddelbart at den største fælles frekvens for 100Hz, 150Hz og 300Hz er 50Hz, hvorfor at Tb = 1/50Hz = 0.02s.

Herunder ses en enkelt periode af g(n) plottet

o2signal.png

Samme funktion som i opgave 1 er benyttet til at finde spektret af g(n).

o2spect.png

Der ses tydelige deltafunktioner i signalets frekvens komponenters, 100Hz, 150Hz og 300Hz, med amplituder der svarer til halvdelen af amplituden af komponenten i signalet. Det bemærkes at for 0Hz komponenten er amplituden lig amplituden af komponenten i signalet. Dette ses også af (5.30) og (5.31) i noterne.

Filteret h(n) er defineret som følger:
$h(n) = 1/N, 1 \leq n \leq N, h(n)=0 \text{ elsewhere}, N=10$
I teori forberedelse opgave 3 er dette filters overføringsfunktion fundet analytisk. Ved at benytte plot_spectrum.m findes dets overføringsfunktion computationelt og de 2 plottes samt signalets spektra G(f)

o2Hspect.png

Der ses ingen forskel mellem den analytiske og den udregnede overføringsfunktion. Det aflæses at den har værdien 1 for 0Hz, 0 for 100Hz, 0.2 for 150Hz og 0 for 300Hz. Signalet forventes altså ændret efter filtrering således at DC komponenten passerer ufiltreret, 100 Hz og 300Hz filtreres fulstændigt fra, og 150 Hz får 0.2 gange så høj amplitude (0.2*0.2 = 0.04). Altså:

o2gfilt.png

Nu foldes 4 perioder af g(n) med h(n), hvilket giver

o2conv.png

Det ses at den første periode, fra 0s til 0.02s er aperiodisk, mens at den tredje periode, fra 0.04s til 0.06s er periodisk og ligner det forventede filtrerede signal, med en DC på 1.2 og en enkelt sinus komponent med 3 perioder på 0.02 sekunder (150Hz) med en amplitude på ca. 0.04 ((1.244-1.161)/2).
Det forventes altså at spekret af den første periode ikke vil ligne spektret der fremkom ved G(f)*H(f), men at spektret af tredje periode vil ligne den i stor udstrækning.

o2first.png

Der ses en tydelig afvigelse fra G(f)*H(f)

o2third.png

Og der ses en tydelig lighed med G(f)*H(f)

MATLAB-koden brugt til denne opgave ses nedenfor:

clear all; close all; clc;
fs = 1000; %Hz - sampling frequency

% The fundamental frequency, fb, is the greatest common divisor of the
% frequncies:
fb = 50; %Hz

% One period of g is therefore:
T = 1/fb;
[n, g] = o2(T); %One period of the signal g
figure;
plot(n./fs,g,'.-');
xlabel('Time (t) [s]'),ylabel('Amplitude [Voltage]')
title('One period of g(n)')

N = 10;
h(1:N) = 1/N;
% The computed transfer function:
[Hfft,f]=plot_spectrum(h,fs,1000);

% The analytic transfer function:
Hcalc = 1/N*sin(pi*f/fs*N)./sin(pi*f/fs).*exp(-j*pi*f/fs*(N+1));
Hampcalc = abs(Hcalc);
figure;
hold on
plot(f,Hfft,'r-',f,Hampcalc,'b--')

[G,f] = plot_spectrum2(g,fs); % The spectrum of the periodic signal 
stem(f,G,'ko')

legend('Computed H(n)','Analytical H(n)','G(f)');
xlabel('Frequency (F) [Hz]'),ylabel('Amplitude [Voltage]');
title('G(F), Analytical H(n) and computed H(n)');

hold off

figure; %The filtered spectrum
[Hfft,f]=plot_spectrum(h,fs,length(g));
stem(f,G.*Hfft,'o')
xlabel('Frequency (F) [Hz]'),ylabel('Amplitude [Voltage]')
title('The spectrum of G(n)*H(n) (filtered)')

figure %Convolution fun
[n,sig] = o2(4*T);
y=conv(h,sig);
plot((1:length(y))./fs,y);
xlabel('Time (t) [s]'),ylabel('Amplitude [Voltage]')
title('The time signal of conv(g(4*T),h(n)) (filtered)')
plot_spectrum2(y(1:0.02.*fs),1000,1);
title('The spectrum of the first period of conv(g(4*T),h(n))')
plot_spectrum2(y(0.04.*fs:0.06.*fs),1000,1);
title('The spectrum of the third period of conv(g(4*T),h(n))')

Opgave 3

Overføringsfunktionen fundet med MATLAB (vha. plot_spectrum.m fra øvelse 1) er sammenlignet med den analytisk udregnede fra teori afsnit 3.2, resultaterne ses i figuren nedenfor. Det ses at forskellen på de to er i størrelsesorden $10^{-15}$, hvilket er omkring MATLABs afrundings nøjagtighed - det er altså højst sandsynligt kun afrundingsfejl der er ansvarlige for forskellen på de to.

opg3HcHfftt.png

Har filtreret g(n) fra spørgsmål 2 på to måder - dels ved en foldning af h(n) og det i MATLAB repræsenterede g(n), og dels ved at gange overføringsfunktionen og spektret for g(n) fundet i opg 2. Amplitude spektret for G, Y (filtreret G) og H opnået ved de to metoder ses nedenfor. På den første figur ses spektret for G og H plottet. Der ud fra forventes det at en filtrering af g resulterer i et y hvor den eneste frekvenskomponent tilbage er én ved 150Hz, med samme amplitude som i g, nemlig 0.2. Denne konklussion drages idet at H netop er 1 ved F=150HZ, men 0 ved F={0,100,300} (pga. både h(n) og g(n) er reelle signaler er deres spektre symmetriske omkring F=0).

Den næste figur viser en filtrering af g ved multiplikation af G og H, og dér ses netop det forventede spektrum af y.

Den sidste figur viser spektret af "y" opnået ved en foldning af h og det i MATLAB repræsenterede g. Dette spektrum svarer ikke til forventningerne, hvilket skyldes at det ikke er det korrekte g der er filtreret, men det i MATLAB repræsenterede g, som er en firkantfunktion multipliceret med g. Derimod var den anden metode korrekt, da spektret af g godt kan findes korrekt i MATLAB ud fra én periode af g.

opg3GH.png
opg3Y.png
opg3yconv.png

MATLAB-koden brugt til denne opgave ses nedenfor:

close all; clear all; clc;

fs=1000; %sampling frequency
% Impulse response of the filter:
h = [0.5, zeros(1,9), -0.5]; 

%% The transfer function:
% Transfer function determined using MATLABs fft:
[Hfft,f]=plot_spectrum(h, fs, 1000);

% Calculated transfer function from theory preperation
Hc = j*sin(10*pi*f/fs).*exp(j*10*pi*f/fs); Hcamp = abs(Hc);

% Comparing the two:
figure;
plot(f,Hfft-Hcamp); 
legend('|Hfft|-|Hc|')
title('Difference between (amplitude spectri of) the  transfer function found with MATLAB (Hfft) and the transfer function calculated analytic (Hc)')
xlabel('Frequency (F) [Hz]'),ylabel('Amplitude [Voltage]')
% The calculated and the matlab transfer functions are equal. 

% The signal g:
fb = 50 %The basic frequncy
T = 1/fb; % The period
g = o2(T); %One period of the signal g

% The expectation:
[G,fg] = plot_spectrum2(g,fs); % The correct spectrum of 
%the periodic signal!  

% The filtration:
Hc = j*sin(10*pi*fg/fs).*exp(j*10*pi*fg/fs);
%Hcamp = abs(Hc);
Yamp = abs(Hc.*G);

% G and H:
figure;
stem(fg,G,'o')
hold on
plot(f,Hfft, 'k')
hold off
legend('|G|','|Hfft|');
title('|G| and |H|')
xlabel('Frequency (F) [Hz]'),ylabel('Amplitude [Voltage]')

% The spectrum of the filtered signal:
figure;
stem(fg,Yamp,'o')
legend('|Y|') 
title('The spectrum of y')
xlabel('Frequency (F) [Hz]'),ylabel('Amplitude [Voltage]')

% Alternative filtration?
y = conv(h,g);
[Y,fy]=plot_spectrum2(y, fs);
figure;
stem(fy,Y)
legend('|"Y"|')
title('The spectrum of y, achived by convolution of MATLB-g and h(n)')
xlabel('Frequency (F) [Hz]'),ylabel('Amplitude [Voltage]')

Opgave 4

Istedet for at aflæse den gennemsnitlige frekvens blev der udarbejdet en automatisk metode til at gøre dette.
Den gennemsnitlige frekvens findes som

(17)
\begin{align} \frac{\sum_{f=-Fg}^{Fg} G(f) \times f}{\sum_{f=-Fg}^{Fg} G(f)} \end{align}

Hvilket giver

o4velo1.png
o4velo2.png

Følgende matlab kode blev brugt til at generere ovenstående:

close all; clear all; clc

% Loading the us-signals:
load velo1.mat
load velo2.mat

Fs = 5e3; %Hz
M = 25000; %Hz 
f0 = 3e6; %Hz
c = 1540; %m/s

% Calculating the spectra:
[VELO1, F1] = plot_spectrum(velo1, Fs, M, 1);
F_m1 = sum(F1.*VELO1)/sum(VELO1); %Hz The mean frequency determined automatically from the spectra
v1 = F_m1*c/(2*f0); %m/s The blood velocity.
x_m1 = round(F_m1);
line([x_m1 x_m1],[0 max(VELO1)],'color','r')
title(['Mean frequency = ' num2str(x_m1) 'Hz. v = ' num2str(v1) 'm/s']);

[VELO2, F2] = plot_spectrum(velo2, Fs, M, 1);
F_m2 = sum(F2.*VELO2)/sum(VELO2); %Hz The mean frequency determined automatically from the spectra
v2 = F_m2*c/(2*f0); %m/s The blood velocity.
x_m2 = round(F_m2);
line([x_m2 x_m2],[0 max(VELO2)],'color','r')
title(['Mean frequency = ' num2str(x_m2) 'Hz. v = ' num2str(v2) 'm/s']);
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License