BME BME 4 4 52 52 Bio Bio medi medi cal Signal cal Signal Processing
Processing Lecture 4 Lecture 4
Fourier analysis and spectral Fourier analysis and spectral estimation
estimation
Lecture 4 Outline
In this lecture, we’ll study the frequency domain analysis methods for signals
Discrete Fourier transform
Classical spectral estimation
Some examples of spectral estimation of signals
Introduction
Knowledge of cyclic or oscillating activity in signals is very important
Fourier transform is the oldest method for analysing the cycles, i.e.
obtaining frequency information
Eg: analysing the cycle of variable stars is the oldest application of analysing cycles
Variable stars are, quite simply, stars that change brightness
An example is given below:
What is the approximate period?
Figure from R.Shiavi, Introduction to applied statistical signal analysis
Discrete Fourier Transform
Since we are dealing with discrete-time signals, we will skip continuous Fourier transform and move into discrete Fourier transform (DFT)
Discrete time and frequency
We have discrete time, similarly, we can have only discrete frequencies
Discrete frequency
Assume N is the number of sampled points
Define the duration of the signal, P=NT, where T is the sampling interval
We have the discrete frequency variable (frequency number), m
The frequency number of m will run from 1 to N
The actual frequency will run from 0 to fs/2 with fd intervals
So, the frequency spacing, fd=1/NT or fd=fs/N
Discrete Fourier Transform (cont)
Some important points
The actual frequency =m*fd=m*fs/N from 0 to fs
Frequency number, m=0,1,2,3,…..N points
Due to symmetry, we normally plot the values from 0 to fs/2 only
Actual frequency vs frequency number
Consider the following example
Time = 2s, N=12 => fs = 6 Hz
Frequency plot is from 0 to 3 Hz
m=0,1...N/2 (the rest of N/2 points rejected due to symmetry in freq domain)
So, m = 0,1,2,3,4,5,6
Actual frequency => (m*fs)/N
So, actual frequency => 0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0
Discrete Fourier Transform (cont)
DFT is used computationally to perform Fourier analysis of discrete-time signals
DFT is defined as
Inverse DFT (IDFT) is defined as
1 0
2
) ( )
(
N
n
N mn j
DFT m T x n e
X
1
0
2
) 1 (
) (
N m
N mn j
e m NT X
n x
Note: i and j will be used to
represent 1
Discrete Fourier Transform (cont)
Using Euler’s relation,
So (why?)
Now, DFT is
) sin(
)
cos(x j x
e jx
) sin(
)
cos(x j x e jx
1
0
sin 2 cos 2
) ( )
(
N
n
DFT N
j mn N
n mn x T
m
X
1
0
1 0
sin 2 ) 2 (
cos ) ( )
(
N n
N n
DFT N
n mn x
N jT n mn
x T
m
X
Discrete Fourier Transform (cont)
As you can see, DFT now consist of a real part and an imaginary for every frequency number m
The magnitude spectra
And phase spectra
where
2
2 Im( ( ))
)) ( Re(
)
(m X m X m
X
)) ( Re(
)) ( arctan Im(
)
( X m
m m X
1
0
cos 2 ) ( )
Re(
N
n N
n mn x
T
m
1
0
sin 2 ) ( )
Im(
N
n N
n mn x
T
m
Normally, we are interested in this for signals
Discrete Fourier Transform using MATLAB
2
2 Im( ( ))
)) ( Re(
)
(m X m X m
X
)) ( Re(
)) ( arctan Im(
)
( X m
m m X
1
0
cos 2 ) ( )
Re(
N
n N
n mn x T
m
1
0
sin 2 ) ( )
Im(
N
n N
n mn x T
m
for m=1:N, sumn=0;
for n=1:N,
sumn=sumn + x(n)*(cos (2*pi*(m-1)*(n-1)/N));
end
dftreal(m)=T*sumn;
end
for m=1:N, sumn=0;
for n=1:N,
sumn=sumn + x(n)*(sin (2*pi*(m-1)*(n-1)/N));
end
dftimag(m)=-T*sumn;
end
for m=1:N,
dftmag(m)=sqrt(dftreal(m)*dftreal(m)+dftimag(m)*dftimag(m));
dftphase(m)=atan(dftimag(m)/dftreal(m));
end
You can also use fft function in MATLAB, do you know what is fft?
A worked example of DFT
The table shows the intensity of ionospheric reflections starting from midnight
Assume fs=24 (i.e. hours in a day)
T=1/24, N=12
Time(hrs) Intensity
0 6
1 20
2 28
3 8
4 1
5 -7
6 20
7 6
8 7
9 -14
10 -19
plot(x);
axis([1 12 -20 30]);
title('Ionospheric data');
xlabel('Time (h) - note it should one less');
ylabel('Intensity');
Figure from R.Shiavi, Introduction to applied statistical signal analysis
MATLAB indexes from 1, so 1 must be reduced from this index to push
Ionospheric DFT
plot(dftreal,'+-');
xlabel('Frequency number (m)');
ylabel('Magnitude');
hold on;
plot(dftimag,'r*-');
axis([1 12 -4 4]);
xlabel('Frequency number (m) - note it should one less');
title('Real and imaginary parts');
Note the symmetry line.
DFT phase is
conjugate symmetry at this line.
Y= a+jb Y*=a-jb
Figure from MATLAB
Figure from R.Shiavi, Introduction to applied statistical signal analysis
Frequency number or actual frequency value??
If the frequency number is m, the actual frequency values are m*fs/N –more on this later
Ionospheric Spectra
For magnitude and phase spectra, the actual frequencies are m*fs/N
1 2 3 4 5 6
-2 -1 0 1 2 3 4
Scale is in m, but note: frequency in cycles/day is 2(m-1)
Magnitude or phase(rad)
Ionospheric Spectra
Figure from MATLAB
Figure from R.Shiavi, Introduction to applied statistical signal analysis
Some points to note
Normally the magnitude spectra is the output that we want from DFT
Magnitude spectra measure the energy of the signal at the specific frequencies
Magnitude spectra runs from f=0 to one point less than fs/2 with N/2 points, i.e. fs/N frequency spacing
In the example, fs=24, N=12, so the actual frequency scale runs from 0,2,4,6,8,10 Hz
But DFT is plotted in frequency numbers for N points, so frequency number runs from 0,1,2,3,4,….,11
Ionospheric Spectra (cont.)
Picket fence effect
One of the problem with DFT is that it provides values of X(m) only at a specific set of frequencies
What if an important value exist at the frequency, f≠m?
In our example, we have magnitude spectra values for freq=0,2,4,6,8,10
What if there was an important cycle with freq=1? We would have missed it!
This problem is known as picket fence effect
Zero padding to increase the DFT resolution
The solution is to increase the signal size by zero padding
This will increase N and hence reduce the frequency spacing
And increase the resolution of DFT
Example, assume we pad 12 zeros to the ionospheric signal and compute the magnitude spectra
So, now N=24
Freq spacing =fs/N=24/24=1
And frequency scale has 12 points and has values for freq=
0,1,2,3,4,5,6,7,8,9,10,11
From the plot, we can see that there are values for freq=1,3,5,
….
Figure from R.Shiavi, Introduction to applied statistical signal analysis
Effect of truncation
When we compute DFT, we assume that the signal is one complete cycle extracted from a portion of a longer signal that is periodic
This is a basic assumption of DFT
That is, we have truncated the signal (though in reality we did not)
If the signal was indeed periodic, it would the same starting and ending points after truncation
But then for our signals, we have discontinuity, so this creates a lot of noise
Technically this problem is known as spectral leakage but we won’t worry about this term
In the ionospheric signal example:
Discontinuity at Truncate Ionospheric signal
Windowing
Now, assume if we have the triangular signal as shown below
Will there be a problem with DFT?
No, because the start and end points are the same
Very often, we consider using windows other than rectangular window to ‘force’ the signal value at the beginning and end to the same value
Normally, this value is zero, though not always
But how we do this?
We use spectral windows like those on the right
Figure from R.Shiavi, Introduction to applied statistical signal analysis
Triangular signal
Windowing (cont.)
The simplest window is the rectangular window
Using this window, the value at beginning (n=0) and end of signal (n=N-1) is set to zero
But this creates discontinuity at second and second last points, i.e. between points n=0 and n=1 and at points n=N-2 and n=N-1
So, if we continue to reduce discontinuity using this window, we land up with a zero signal!
To solve this, we can use smoothed windows like Triangular, Blackman, Hamming, etc.
The new windowed signal for every point n, y[n]=x[n]*w[n]
Some window functions:
Blackman:
Hamming:
Triangular:
For these functions, n is the window size (i.e. N)
These window equations will not be tested
Windowing (examples)
The following show the magnitude spectra of an EEG signal:
0 20 40 60 80 100 120 140
0 10 20 30 40 50 60
0 20 40 60 80 100 120 140
0 10 20 30 40 50 60
0 20 40 60 80 100 120 140
0 5 10 15 20 25 30
Original, no windowing
Rectangular windowing
Hamming windowing
Windowing gives a smoother DFT
The choice of windowing is too complex for this course, it involves analysis of ratios of magnitudes main lobes and side lobes
So, we’ll skip them and choose any window (other than rectangular)
Do you know the sampling frequency of the EEG from the plot?
Detrending
Figures from R.Shiavi, Introduction to applied statistical signal analysis
Another aspect similar to windowing is the removal of polynomial time trends and average value
The average value, which has energy around zero frequency can be removed simply by setting the mean to zero
The polynomial time trends have energy at very low frequencies
We’ll study some high pass filtering methods in the next lecture to remove this polynomial trend
This is known as detrending. Look at the example below:
Before detrending After detrending
DFT applied to some artificial signals
Consider the sinusoidal wave with frequency=10 and fs=200
N=200;
X=sin(2*pi*(1:N)*10/200);
plot(X);
The signal and magnitude spectra are shown below
The same signal but length N=20
0 50 100 150 200
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
0 20 40 60 80 100
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Do you need windowing for this signal?
What has happened – what is the actual frequency now?
Since fd=fs/N=10,
MATLAB indexes from 1, so 1 must be reduced from this index to push it back to 0
Freq(Hz) =(m-1)*fs/N =(2-1)*10 Hz Theoretically, is there anything wrong with
computing Fourier spectra for this signal this way if fs is known?
DFT applied to some artificial signals (cont.)
Consider combination of sinusoidal waves with frequency=8 and 25 and fs=200
N=25;fs=200;f1=8;f2=25;
X=sin(2*pi*(1:N)*f1/fs)+sin(2*pi*(1:N)*f2/fs);
plot(X);
The signal and magnitude spectra are shown below
Theoretically, the minimum N for this signal must N=25.
Why?
The same signal but length N=200
The magnitude spectra does show two peaks but there is a problem – what is it?
0 5 10 15 20 25
-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
0 2 4 6 8 10 12
0 0.2 0.4 0.6 0.8 1 1.2 1.4
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-1 -0.5 0 0.5 1 1.5
2 - Peak at 9 and 26 =>Spectra
at 8 and 25 Hz
- Both have same magnitude In reality, DFT works properly only if there are multiple cycles and not just one cycle And there must enough
MATLAB indexes from 1, so 1 must be reduced from this index to push it back to 0
Remember: Freq (Hz)=(m-1)*fs/N
DFT- a worked example
Consider the following frequency spectra (the frequency axis should be up to 500 but is shown partially only)
Given N=1000 and fs=200, find out the values for f1, f2, A1, A2 for
MATLAB code : X=A1*sin(2*pi*(1:N)*f1/fs)+A2*sin(2*pi*(1:N)*f2/fs);
Answer
Freq (Hz)=(m-1)*fs/N
A1=1.0
A2=0.5
f1=10
f2=18
The ratio of A1/A2 has to be 2 but the exact values could be different
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
The peaks are at 51 and 91
There are 3 common PSD estimation using classical approaches
Periodogram
Welch method
Using autocorrelation function (but will not be studied here)
There is a recent advanced parametric method on PSD
estimation based on autoregressive modelling but this will not be studied here.
Power spectral density estimation, classical approaches
What is Power Spectral Density?
PSD is similar to the square of magnitude spectra
Actually, the only difference is the scale factor
So, PSD=power per unit frequency (in dB)
If we compute magnitude spectra using DFT, convert it to PSD, then we have periodogram
Sometimes known as Blackman-Tukey method
You can use the codes below OR use psd(x) or periodogram(x) functions
0 20 40 60 80 100
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1 N=200;fs=100;f1=5;
X=sin(2*pi*(1:N)*f1/fs);
plot(X);
Y=fft(X,N);
Pyy = sqrt((Y.* conj(Y))) /(0.5*N);
plot(Pyy(1:N/2));
PSD=10*log10(Pyy.*Pyy);
plot(PSD(1:N/2))
0 20 40 60 80 100
-350 -300 -250 -200 -150 -100 -50 0
Computing PSD using Welch method
An improvement to periodogram
Normally gives a smoother PSD
The method
Segmentation with overlap
The sample data is divided into K segments each containing M sample points
Initial data points in each
segment are separated by D (i.e.
non-overlap length is D)
Segments overlap with M-D samples and K=(N-M)/D +1 segments are formed
The percent overlap is 100*(M- D)/M
Apply the conventional periodogram (i.e. DFT) to compute PSD
In MATLAB, use pwelch(x)
function where the percentage overlap is 50% and Hamming window is used
Computing PSD using Welch method (example)
0 0.2 0.4 0.6 0.8 1
-10 -5 0 5 10 15 20 25 30 35
Frequency
Power Spectrum Magnitude (dB)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
-15 -10 -5 0 5 10 15 20 25 30
Normalized Frequency ( rad/sample)
Power/frequency (dB/rad/sample)
Power Spectral Density Estimate via W elch
Assume we have an EEG signal with N=125
The first plot is obtained by using psd, while the second is obtained by using pwelch with each segment, M=50
D=25 (since overlap percentage is 50%)
So, overlap of M-D=25 points and K=(125-50)/25+1=4
From the plots, pwelch is more smoother and very often more accurate than psd/periodogram
PSD frequency scale is normalised and in radian instead of cycles
0 to fs 0 to 2 in (Hz rad/sample) So, 0 to fs/2 0 to
Homework 1 - DFT (worked example)
Compute the DFT for the ionospheric data and see if you get the values as in the table
Next compute the magnitude and phase spectra and see if you get the values in the plots as below
Homework 2
For the following noise corrupted signal, window the signal using triangular window, then compute the magnitude spectra (using MATLAB)
Can you get similar plots as below?
Use the codes below (fft is used)
N=100;fs=200;f1=35;
X=sin(2*pi*(1:N)*f1/fs)+rand(1,N);
plot(X);
0 20 40 60 80 100
-1 -0.5 0 0.5 1 1.5 2
0 20 40 60 80 100
-1 -0.5 0 0.5 1 1.5 2
W=triang(N);
X=X'.*W;
plot(X);
Y=fft(X,N);
Pyy = sqrt((Y.* conj(Y))) /(0.5*N);
plot(Pyy(1:N/2));
0 10 20 30 40 50
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Study guide (Lecture 4)
From this week’s lecture, you should know how to compute
Discrete Fourier Transform
PSD estimation using classical methods (periodogram, Welch)
You will need to remember all the equations in this lecture
Note: the use of MATLAB codes will never be tested in the exam. For example, you will never be tested on what is fft(x), what is psd(x), etc. or how to use them.
End of lecture 4
