Digital Representation of Audio Information - University of Kentucky

pancakesbootΤεχνίτη Νοημοσύνη και Ρομποτική

24 Νοε 2013 (πριν από 3 χρόνια και 11 μήνες)

74 εμφανίσεις

EE513

Audio Signals and Systems

Digital Signal Processing (Analysis)

Kevin D. Donohue

Electrical and Computer Engineering

University of Kentucky

Spectra


The spectrum of the impulse response indicates the following
filter/system properties:


resonance and frequency sensitivities through the magnitude spectrum


delay and dispersive properties through the phase spectrum.



The discrete Fourier transform (DFT) computes complex
samples of the spectrum.


For deterministic signals both the phase and magnitude are important
for characterizing the signal or response.


For stationary noise processes the square of the DFT magnitudes can
be averaged from independent segments to create a power spectral
density.

Discrete Fourier Transform


In order to take a DFT of a signal, a finite window (time interval) must be
extracted from the original sampled time signal:










The process can be repeated with overlapping windows covering entire signal.


If the signal is deterministic and changes over time, the DFTs for small consecutive intervals
must be displayed to show the changes in spectrum over time (Spectrogram).


If the signal is random with statistics that do not change over time, the DFT magnitudes can
be averaged (Power Spectral Density (PSD)).

3.2
3.4
3.6
3.8
4
4.2
4.4
4.6
4.8
5
5.2
x 10
4
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
Samples
4.34
4.345
4.35
4.355
4.36
4.365
4.37
4.375
4.38
4.385
4.39
x 10
4
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
Samples
-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
0
200
400
600
800
1000
Hz (Normalized Fs = 1)
DFT Magnitude Squared
-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
-4
-3
-2
-1
0
1
2
3
4
Hz (Normalized Fs = 1)
DFT Phase in Radians
DFT Phase

DFT and Inverse DFT



DFT points are evaluated at discrete frequency samples
(frequency sampling). The DFT for an N point sequence is
defined as:






This inverse DFT is a similar computation to the DFT
except the complex conjugate of the kernel is used:



N
n
N
n
jk
k
X
N
n
x
N
k












0
for

2
exp
)
(
ˆ
1
1
0



N
k
N
k
jn
n
x
k
X
N
n













0
for

2
exp
)
(
ˆ
1
0

DFT with
W

notation


Let
W
N

equal the complex sinusoid with the lowest nonzero
frequency (
k

= 1):



N
k
W
n
x
k
X
N
n
kn
N






0
for

)
(
ˆ
1
0








N
j
W
N

2
exp

W
N

is effectively a fundamental frequency in a harmonic
expansion, since the frequency domain is sampled in equal
increments. DFT and its inverse can now be written as:



N
n
W
k
X
N
n
x
N
k
nk
N







0
for

)
(
ˆ
1
1
0
Computation Examples


Compute the DFT for the 4 point and 8 point
signals given below:





0

,

1,1

,
0
3

,
2

,
1

,
0

)
(
)
(
)
(
)
(
x
x
x
x




0

0,

0,

,
0

0,

,

1,1

,
0
7


,
1

,
0

)
(
)
(
)
(
x
x
x

What was the effect of extra zeros padded on
the end of the signal in the second example?

Computation Examples


Compute the DFT for
the 4 point



0

,

1,1

,
0


4
0
for

)
(
ˆ
3
0
4





k
W
n
x
k
X
n
kn


2
0
1
1
0

)
(
ˆ
,
0
0
4
0
4
0
4
0
4
3
0
0
4









W
W
W
W
W
n
x
k
X
k
n



135
2
1
0
1
1
0

)
(
ˆ
,
1
3
4
2
4
1
4
0
4
3
0
1
4














j
W
W
W
W
W
n
x
k
X
k
n
n


0
0
1
1
0

)
(
ˆ
,
2
2
4
0
4
2
4
0
4
3
0
2
4









W
W
W
W
W
n
x
k
X
k
n
n



135
2
1
0
1
1
0

)
(
ˆ
,
3
1
4
2
4
3
4
0
4
3
0
3
4













j
W
W
W
W
W
n
x
k
X
k
n
n
30

210

60

240

90

270

120

300

150

330

180

0

W

3

4


= j

W

2

4


=
-
1

W

0

4

= 1

W

2

4

=
-
j

RE

IM

Computation Examples


Compute the DFT for the
8 point



0

0,

0,

0,

0,

,

1,1

,
0


8
0
for

)
(
ˆ
7
0
4





k
W
n
x
k
X
n
kn


2
0
0
1
1
0

)
(
ˆ
,
0
0
8
0
8
0
8
0
8
7
0
0
8










W
W
W
W
W
n
x
k
X
k
n



5
.
67
85
.
1
71
.
1
71
.
0
0
1
1
0

)
(
ˆ
,
1
3
8
2
8
1
8
0
8
7
0
1
8














j
W
W
W
W
W
n
x
k
X
k
n
n


0
0
0
1
1
0

)
(
ˆ
,
4
4
8
0
8
4
8
0
8
7
0
4
8










W
W
W
W
W
n
x
k
X
k
n
n



135
2
1
.
0
0
1
1
0

)
(
ˆ
,
2
6
8
4
8
2
8
0
8
7
0
2
8















j
W
W
W
W
W
n
x
k
X
k
n
n



5
.
157
76
.
0
29
.
0
71
.
0
0
1
1
0

)
(
ˆ
,
3
1
8
6
8
3
8
0
8
7
0
3
8














j
W
W
W
W
W
n
x
k
X
k
n
n



5
.
157
76
.
0
29
.
0
71
.
0
0
1
1
0

)
(
ˆ
,
5
7
8
2
8
5
8
0
8
7
0
5
8















j
W
W
W
W
W
n
x
k
X
k
n
n



135
2
1
0
0
1
1
0

)
(
ˆ
,
6
2
8
4
8
6
8
0
8
7
0
6
8














j
W
W
W
W
W
n
x
k
X
k
n
n



5
.
67
85
.
1
71
.
1
71
.
0
0
1
1
0

)
(
ˆ
,
7
5
8
6
8
7
8
0
8
7
0
7
8













j
W
W
W
W
W
n
x
k
X
k
n
n
30

210

60

240

90

270

120

300

150

330

180

0

W

0

8

= 1

W

7

8

=

W

4

8


=
-
1

W

3

8


=

W

1

8


=

W

2

8

=
-
j

W

6

8


= j

W

5

8


=



j



1
2
2


j

1
2
2


j


1
2
2


j


1
2
2
Computation Examples


What symmetries exist in the W
N

factor in the DFT
computation that suggest a more efficient computation than
a direct multiply and summation? The implementation that
takes advantage of these symmetries is called a Fast Fourier
Transform (FFT).


If the sequence is a power of 2, the FFT is most efficient.
A direct implementation of the DFT requires N
2

complex
multiplications. An FFT requires Nlog
2
N complex
multiplications.


Compute the ratio between complex multiplications for the
FFT and DFT for a 128 point sequence. Compute the
same ratio for a 2048 point sequence.


See
http://www.cmlab.csie.ntu.edu.tw/cml/dsp/training/coding/transform/fft.html

For more details on the FFT derivation

Homework(1)


Derive each of the symmetry relationships
below for the
W
N

factor:



1

and

,
any
for

*
)
(





N
k
n
W
W
k
n
N
N
nk
N
even


and

,
1
for

2
N
N
W
W
k
N
N
k
N










even


and

1

,
,
any
for

2
2
N
N
k
n
W
W
nk
N
nk
N







1

,
any
for

1
0
*







N
k
N
W
W
N
n
nk
N
nk
N


1
for




N
W
W
k
N
N
k
N
Frequency Sampling, FS to DFT


Time sample the following FS equations to show it
equivalence to the DFT.






Sample time axis with
T
s

where
NT
s

=
T

o
k
o
f
T
T
t
x
t
x
t
kf
j
k
X
t
x
1
),
(
)
(
for

)
2
exp(
]
[
ˆ
)
(












T
o
dt
t
kf
j
t
x
T
k
X
)
2
exp(
)
(
1
]
[
ˆ

DFT Harmonics (N=4)

0
1
2
3
4
0
0.2
0.4
0.6
0.8
1
Samples
x(n) with kernel for k= 0
0
1
2
3
4
-1
-0.5
0
0.5
1
Samples
x(n) with kernel for k= 1
0
1
2
3
4
-1
-0.5
0
0.5
1
Samples
x(n) with kernel for k= 2
0
1
2
3
4
-1
-0.5
0
0.5
1
Samples
x(n) with kernel for k= 3
DFT Harmonics (N=8)

0
2
4
6
8
0
0.2
0.4
0.6
0.8
1
Samples
x(n) with kernel for k= 0
0
2
4
6
8
-1
-0.5
0
0.5
1
Samples
x(n) with kernel for k= 1
0
2
4
6
8
-1
-0.5
0
0.5
1
Samples
x(n) with kernel for k= 2
0
2
4
6
8
-1
-0.5
0
0.5
1
Samples
x(n) with kernel for k= 3
0
2
4
6
8
-1
-0.5
0
0.5
1
Samples
x(n) with kernel for k= 4
0
2
4
6
8
-1
-0.5
0
0.5
1
Samples
x(n) with kernel for k= 5
0
2
4
6
8
-1
-0.5
0
0.5
1
Samples
x(n) with kernel for k= 6
0
2
4
6
8
-1
-0.5
0
0.5
1
Samples
x(n) with kernel for k= 7
FFT Spectral Estimates


Distortion and error result from the truncated
and frequency sampled signal because:



Lost energy from truncated time portion will
distribute error over the signal.


Truncating signal at the endpoints can introduce a
sharp transition that may not be part of the signal.


Sampling signal in frequency produces periodicity
in time, which results in
circular

time
-
domain
convolution from multiplying spectra in frequency
domain.

Frequency Sampling Effect


The DFT can be applied to a finite duration signal. If signal support is
infinite, part of it must be truncated. The DFT definition is applied to the
finite number of points to obtain a frequency sampled spectrum, which
results in a periodic extension of the signal in the time window.

0
5
10
15
20
25
30
35
40
45
50
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Sampled signal for DFT analysis
Time in Seconds
Amplitude
Truncate Signal 30 Seconds
0
50
100
150
200
250
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Sampled signal for DFT analysis
Time in Sample Index
Amplitude
Truncate Signal 150 Samples
Frequency Sampling Effect


The DFT representation implies the time domain signal was
periodic (sampling in time produces periodic spectra, likewise
sampling in frequency produces periodic time segments) as
would be the case for Fourier series.

-150
-100
-50
0
50
100
150
200
250
300
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Periodic Extension from Frequency Sampling
Time in Sample Index
Amplitude
Original
150 point signal
Circular Convolution


Multiplication in frequency is convolution in time. For
discrete signals, this can be performed with the fft
command in Matlab on the filter coefficients
h
(
n
) and input
signal
x
(
n
), then the ifft command on their frequency
domain product:







0
m
m
n
h
m
x
n
h
n
x
n
y
)
(
)
(
)
(
*
)
(
)
(
)
(
ˆ
)
(
ˆ
)
(
ˆ
k
H
k
X
k
Y







1
0
)
(
ˆ
N
n
kn
N
W
n
x
k
X






1
0
)
(
ˆ
N
n
kn
N
W
n
h
k
H







1
0
)
(
ˆ
1
N
k
nk
N
W
k
Y
N
n
y
Time

Domain

Frequency

Domain

Time

Domain

Linear Convolution Example


Linear
convolution
implemented
in the time
domain.
Note the
convolved
signal length
is greater
than either of
the originals.
Why?

20
40
60
80
100
120
140
-0.1
0
0.1
0.2
0.3
Filter Impulse Response 128 Points
20
40
60
80
100
120
140
-4
-2
0
2
Signal 128 Points
20
40
60
80
100
120
140
-1.5
-1
-0.5
0
0.5
Linear Convolution in Time Domain 148 Points
Circular Convolution Example


Convolution
implemented
in the
frequency
domain.
Note the
convolved
signal length
is the same as
the originals.
Why?


Why is there
an apparent
artifact at the
beginning of
the signal?

20
40
60
80
100
120
140
-0.1
0
0.1
0.2
0.3
Filter Impulse Response 128 Points
20
40
60
80
100
120
140
-4
-2
0
2
Signal 128 Points
20
40
60
80
100
120
140
-1.5
-1
-0.5
0
0.5
Circular Convolution Implemented in Frequency Domain 128 Points
Circular Convolution Example


Convolution
implemented
in the
frequency
domain with
zero padding.
Why is the
apparent
artifact from
the nonzero
pad example
no longer
present?

50
100
150
200
250
-0.1
0
0.1
0.2
0.3
Filter Impulse Response 256 Points by Zero Padding
50
100
150
200
250
-4
-2
0
2
Signal 256 Points by Zero Padding
50
100
150
200
250
-1.5
-1
-0.5
0
0.5
Circular Convolution Implemented in Frequency Domain 256 Points
Homework(2)


Use the FFT to filter a 32 point square pulse
s
[
n
]:




(index
n

starts a 0)


where the FIR filter has impulse response
h
[
n
]:




a)
Use the fft command to take signals into frequency domain without
padding with zeros. Take inverse fft to obtain time domain signal. Plot the
filtered signals and explain what you observe.

b)
Use the fft command and pad with zeros to double the signal lengths and
repeat part (a)





elsewhere

0
16
29
for

1
]
[






n
n
s




12
32
for

0
0
12
for

6
.
]
[









n
n
n
h
n
FFT for Signal Analysis


The FFT can be applied directly to finite duration
energy signals to examine the spectral content of
the signal.


For Matlab’s FFT:



FFT(X) is the discrete Fourier transform (DFT) of
vector X. For matrices, the FFT operation is
applied to each column.


FFT(X,N) is the N
-
point FFT, padded with zeros
if X has less than N points and truncated if it has
more.

FFT of a Rayleigh Ring


Write a Matlab script to Generate a tone
modulated by a Rayleigh envelope. Plot the
signal, its spectra, and play the sound.

% Simulate a Rayleigh envelope ring, truncate it, take its FFT and examine

% its spectral content


fr = 250; % Frequency of the ring

sigma = .05; % effective duration TIMES 2 of ring

fs = 8000; % Sampling frequency

dur = .25; % signal duration in seconds

nfft = 4096;

t = [0:round(dur*fs)
-
1]/fs; % Create time axis

ring = (t/sigma).*exp((
-
t.^2)/(sigma^2)).*sin(2*pi*t*fr); % Generate Ring signal


figure(1)

plot(t, ring) % Plot ring

title('Rayleigh Envelope Ring Signal')

xlabel('Seconds')

ylabel('Amplitude')

soundsc(ring,fs) % Play sound

FFT of a Rayleigh Ring

spec0 = fft(ring); % No Zero pad

spec1 = fft(ring,nfft); % Zero pad (or truncate) to NFFT points

spec2 = fft(ring,2*nfft); % Zero pad to 2*NFFT


% Create frequency axes for each of the spectra

faxis0 = fs*[0:length(spec0)
-
1]/length(spec0);

faxis1 = fs*[0:length(spec1)
-
1]/length(spec1);

faxis2 = fs*[0:length(spec2)
-
1]/length(spec2);

% Plot spectrum with zero padded version on the same graph and compare

figure(2) % Magnitude

plot(faxis0, abs(spec0),'k',faxis1, abs(spec1),'r
--
',faxis2, abs(spec2),'g.')

legend([int2str(length(t)) ' point FFT'], [int2str(nfft) ' point FFT'], [int2str(2*nfft) ' point FFT'])

title('Magnitude ffts with zero padding')

xlabel('Hertz')

ylabel('Magnitude')


figure(3) % Phase

plot(faxis0, 180*phase(spec0)/pi,'k',faxis1, 180*phase(spec1)/pi,'r
--
',faxis2, 180*phase(spec2)/pi,'g.')

legend([int2str(length(t)) ' point FFT'], [int2str(nfft) ' point FFT'], [int2str(2*nfft) ' point FFT'])

title('Phase ffts with zero padding')

xlabel('Hertz')

ylabel('Degrees')


FFT of a Rayleigh Ring

210
220
230
240
250
260
270
280
290
0
10
20
30
40
50
60
70
80
90
100
ffts with zero padding
Hertz
Amplitude
2000 point FFT
4096 point FFT
8192 point FFT
0
0.05
0.1
0.15
0.2
0.25
-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
0.2
0.3
0.4
0.5
Rayleigh Envelope Ring Signal
Seconds
Amplitude
200
210
220
230
240
250
260
270
280
290
300
-250
-200
-150
-100
-50
0
50
Phase ffts with zero padding
Hertz
Degrees
2000 point FFT
4096 point FFT
8192 point FFT
(Zero Padding Effects)

Homework(3)


Create a signal consisting of 2 sine waves with amplitude 1
and sampled at 4000 Hz. Set one with frequency 250 Hz
and the other with 254 Hz.


a) Take the DFT using window length of 0.05, 0.25 and
0.5 seconds. Describe what you see and make
generalizations about the impact of signal length on
frequency resolution.


b) Repeat part (a) using zeros padding so that each signal
length is effectively 1 seconds. Describe what you see and
make generalization about the impact zero padding on
frequency resolution.

FFT of Extended Non
-
Stationary Signals


When the signal is long with changing dynamics,
as would be the case for a speech or music signal,
the FFT must be applied to short segments of
these signals.


To reduce truncation effects, a tapering is used to
bring the ends of the segment close to zero a the
points of truncation.


To ensure that information is not lost in the
truncation and tapering, an overlap between signal
segments is used.

Impact of Windowing


The windowed DFT is the
product of the original
signal and the windowing
function:



where
x
[
n
] is a signal of
infinite extend,
n

is the
sample index,
g
[
n
] is a
finite widow.



]
[
]
[
]
;
[
m
n
g
n
x
m
n
x
w


0
100
200
300
400
500
-1.5
-1
-0.5
0
0.5
1
1.5
Samples
Amplitude
Hamming Window 2
Hamming Window 1
0
100
200
300
400
500
-3
-2
-1
0
1
2
3
samples
amplitude
Impact of Windowing


In the frequency domain this
becomes convolution with the
window DFT:

)
2
exp(
]
[
ˆ
]
[
ˆ
]
;
[
ˆ
1
0
km
j
l
k
G
l
X
m
k
X
N
l
w






0.02
0.04
0.06
0.08
0.1
0
0.05
0.1
0.15
0.2
0.25
0.3
Hertz
Magnitude
Not windowed (very long segment)
Boxcar window
Hamming window
0
100
200
300
400
500
-4
-2
0
2
4
samples
amplitude
Hamming window
0
100
200
300
400
500
-4
-2
0
2
4
Boxcar window
0
0.1
0.2
0.3
0.4
0.5
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
Hertz
Magnitude
Hamming window
Not windowed (very long segment)
Boxcar window
Spectrogram for Signal Analysis


For a signal that deterministically changes spectral properties over
time, the DFT over local time windows must be computed and
indexed. This function is referred to as the
spectrogram:




where
m

is the time index, T is the sample increment to the next
window position, which effectively slides the
N

point analysis window
over consecutive (often overlapping) segments of the signal
x
[
n
]. The
DFTs are plotted along the y
-
axis, while the x
-
axis is time (usually
corresponding to the short
-
time window centers. The result is a
function that shows how the spectrum changes over time.


The following parameters must be set properly for a usable analysis:


window length N, window type, window increment (overlap between
consecutive windows), number of FFT points.







1
0
]
[
]
[
]
;
[
ˆ
N
n
nk
N
N
W
n
g
mT
n
x
m
k
X
Spectrogram for Signal Analysis


The Matlab
spectrogram

command is appropriate for analyzing long
deterministic or non
-
stationary signals.



For Matlab’s SPECTROGRAM:


>> S = SPECTROGRAM(X,WINDOW,NOVERLAP,NFFT,Fs)


calculates the spectrogram for the signal sampled at
Fs

in vector
X
.
SPECTROGRAM splits the signal into overlapping segments,
windows each with the
WINDOW

vector and forms the columns of
S

with their zero
-
padded, length
NFFT

discrete Fourier transforms.
Thus each column of
S

contains an estimate of the short
-
term, time
-
localized frequency content of the signal
X

(positive frequencies only).
Time increases linearly across the columns of
S
, from left to right.

Spectrogram for Signal Analysis

If plotting, it is convenient to use the optional output arguments:




[S,F,T] = SPECTROGRAM(A, WINDOW, NOVERLAP, NFFT, Fs)
returns a column of frequencies
F

for each time in
T

at which the
spectrum was computed.
F

has length equal to the number of rows of
S
, and
T

has length equal to the number of columns. If you leave
Fs

unspecified, SPECTROGRAM assumes a default of 2 Hz. Optional
output parameters
F

and
T

are the frequency and time axes useful for
plotting the spectrum:


>> imagesc(T, F, abs(S))


Spectrogram of Frequency Sweep


Write a Matlab script to Generate a frequency sweep from
20 Hz to 1.9 kHz, where Fs is 4 kHz. Then plot the
magnitude of the spectrogram and play the sound.

% Simulate a frequency sweep signal and take its spectrogram and examine

% its spectral content over time


% Signal Parameters

flow = 20; % Starting frequency of sweep in Hertz

fend = 1900; % Ending frequency of sweep in Hertz

dt = 4; % Time duration of sweep in seconds

fs = 4000; % Sampling frequency


% Spectrogram Parameters

wlen = 128; % Length of actual point extracted from signal segment

nfft = 1024; % number of FFT point (zero padding)

olap = floor(wlen/2); % Points of overlap between segments

wn = hamming(wlen); % Create tapering window (also try boxcar
-

square window)

Spectrogram of Frequency Sweep

% Generate signal

t = [0:round(dt*fs)
-
1]/fs; % Create time axis

fsw = flow + ((fend
-
flow)/2)*[0:length(t)
-
1]/length(t); % Create Frequency ramp
-

Why divid by 2?

swp = sin(2*pi*t.*fsw); % Generate sweep signal

soundsc(swp,fs) % Play sound


% Create Spectrogram

[b,faxis,taxis] = spectrogram(swp,wn,oplap,nfft,fs);

% Plot over time and frequency

figure(1); imagesc(taxis, faxis, abs(b)) % Plot spectrogram

axis('xy') % Flip y axis to put zero Hz on bottom

colorbar % Include colorbar to determine color coded magnitudes on graph


title(['Frequency Sweep from ' num2str(flow) ' Hz to ' num2str(fend) ' Hz']); xlabel('Seconds'); ylabel('Hz')


% Just plot at a single time instant


figure(2); tindex = find(taxis >= dt/2); %Find index of halfway point over the time axis


tindex = tindex(1);


plot(faxis,abs(b(:,tindex)))


title(['Single column from Spectrogram at ' num2str(taxis(tindex)) ' seconds']); xlabel('Hz'); ylabel('Magnitude')

Spectrogram of Frequency Sweep

Seconds
Hz
Frequency Sweep from 20 Hz to 1900 Hz
0
0.5
1
1.5
2
2.5
3
3.5
0
200
400
600
800
1000
1200
1400
1600
1800
2000
5
10
15
20
25
30
35
0
200
400
600
800
1000
1200
1400
1600
1800
2000
0
5
10
15
20
25
30
35
Single column from Spectrogram at 2 seconds
Hz
Magnitude
Homework(4)


Use the scale program created in a previous homework problem to
generate a scale for 2 octaves starting at 256 Hz with tones of
amplitude 1 and duration 0.25 seconds. Use a sampling rate of 4kHz.

A) Compute and plot the spectrogram
magnitude in dB
, labeling all axes
correctly. Use the parameters you think best for estimating the
spectrogram for this signal. Comment on the frequencies of the tones
generated with your program and those identified through the
spectrogram.

B) Quadruple the number of FFT points from that used in (A). Compute
and plot the spectrogram, and explain difference observe from part
(A).

C) Increase the window length by a factor of 4 and set the number of FFT
points to twice the amount of the window length. Compute and plot
the spectrogram, and explain the observed changes.


Random Signal Analysis


If a process is changing over time, while the statistics of
the process do not change over time, the process is referred
to as a stationary process. In this case averaging and
moment computations result in useful characterization.



If the signal’s (or process’) probability density function
never changes over time, then it is referred to as strict
-
sense stationary (SSS), if the first and second moments
remain constant over time, the process is referred to as
wide
-
sense stationary (WSS).

Power Spectral Density


The Power Spectral Density (PSD) is the average power in
a signal as a function of frequency. This can be estimated
by averaging DFT magnitudes squared over segments of
signals from the same noise process.





where E[] is the expected value,
S
(
f
) is the data spectrum
plus noise (i.e. modeled as a random variable) and is
the DFT or short
-
time FT estimated from independent
signal segments. is the true or underlying spectrum.







N
n
i
f
S
N
f
S
f
P
1
2
2
)
(
ˆ
1
)
(
E
)
(
)
(
ˆ
f
S
i
)
(
f
P
PSD Example Colored Noise


A popular approach to spectral estimation is to is to
window a long segment of data into a sequences of shorter
windows (similar to the spectrogram) and average all the
FFT magnitudes together. This is referred to as Welch's or
the hopping window method. The Matlab command
pwelch() implements this method.


Example:

Filter white noise through a band
-
pass filter and
compute the PSD. Compare PSD to the transfer function
magnitude of the filter. Assume a sampling rate of 8kHz
and a 6
th

order Butterworth filter with passband from 500
to 1500 kHz.

PSD Example Colored Noise

fs = 8000; % Sampling Frequency

dur = 20; % Sound duration in seconds

ord = 6;

plen = 32; % PSD segment length

% Bandlimit on the filter

f1 = [500]; % lower bandlimit

f2 = [1500]; % corresponding upper bandlimit

% colors for the plots

col = ['g', 'r', 'b', 'k', 'c', 'b'];

no = randn(1,round(fs*dur)); % Generate noise signal

[b,a] = butter(ord,2*[f1 f2]/fs); % Generate filter

% perform filter operation to get colored noise

cno = filter(b,a,no);

% Compute PSD of noise

[p, fax] = pwelch(cno,hamming(plen),fix(plen/2),2*plen, fs);

figure(1); lh = plot(fax,abs(fs*p/2),col(1)) % Plot PSD

set(lh,'LineWidth',2) % Make line thicker

hold on

% Find filter transfer function

[h,fq] = freqz(b,a,2*plen,fs);

plot(fq,abs(h).^2,col(2))

hold off

% Label figure

xlabel('Hertz', 'Fontsize', 14); ylabel('PSD', 'Fontsize', 14)

0
500
1000
1500
2000
2500
3000
3500
4000
0
0.2
0.4
0.6
0.8
1
1.2
1.4
Hertz
PSD
Window Length = 32


Spectral Estimate
Filter Magnitude Response
0
500
1000
1500
2000
2500
3000
3500
4000
0
0.2
0.4
0.6
0.8
1
1.2
1.4
Hertz
PSD
Window Length = 128


Spectral Estimate
Filter Magnitude Response
0
500
1000
1500
2000
2500
3000
3500
4000
0
0.2
0.4
0.6
0.8
1
1.2
1.4
Hertz
PSD
Window Length 512


Spectral Estimate
Filter Magnitude Response
Autocorrelation Function


As a measure of how quickly a signal changes on
average, or for the detection of periodicities, the
autocorrelation can be used. It directly measures the
correlation between samples separate in time by a
particular interval or
lag
.



Define the autocorrelation of a sequence as:






1
0
N
n
n
x
k
n
x
k
r
)
(
)
(
)
(
N
k
N
N-
n
n
x







and

1
n

and

0
for

0
)
(

where
Autocorrelation Example

Repeat previous example for signal generation and
compute the autocorrelation function. Show
autocorrelation for white noise (filter input) and colored
noise (filter output).

% Compute auto correlation of input sequence

mxlag = fs*.02; % Only compute lags up to .1 seconds

[acwno, lagwno] = xcorr(no, mxlag, 'coef');

% Plot lags

figure(2); plot(1000*lagwno/(fs),acwno)

xlabel('milliseconds','Fontsize', 14); ylabel('Correlation coefficient','Fontsize', 14)



% Compute auto correlation of output sequence

mxlag = fs*.02; % Only compute lags up to .1 seconds

[accno, lagcno] = xcorr(cno, mxlag, 'coef');

% Plot lags

figure(3); plot(1000*lagcno/(fs),accno)

xlabel('milliseconds','Fontsize', 14); ylabel('Correlation coefficient','Fontsize', 14)


Autocorrelation Example

Repeat previous example for signal generation and
compute the autocorrelation function. Show
autocorrelation for white noise (filter input) and colored
noise (filter output).

White noise AC

Colored noise AC

-20
-15
-10
-5
0
5
10
15
20
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
milliseconds
Correlation coefficient
-20
-15
-10
-5
0
5
10
15
20
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
milliseconds
Correlation coefficient
Autocorrelation and PSD

The AC and PSD are DFT pairs. The DFT of the AC is
shown below:

DFT of AC

Colored noise AC

-20
-15
-10
-5
0
5
10
15
20
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
milliseconds
Correlation coefficient
0
500
1000
1500
2000
2500
3000
3500
4000
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
Hertz
Magnitude
Homework(5)


Generate 3 seconds of white noise at a sampling rate of 44.1 kHz. Use
the following coefficients in an IIR filter to filter the white noise and
generate pink noise (
bnum

are numerator and
aden

are denominator
coefficients).

bnum = [ 0.04957526213389,
-
0.06305581334498, 0.01483220320740 ]

aden = [ 1.00000000000000,
-
1.80116083982126, 0.80257737639225]

Coefficients from
http://music.columbia.edu/pipermail/music
-
dsp/2001
-
November/046099.html

A)

Estimate the AC and PSD of the Pink noise and present their plots.
Comment on the observed differences from white noise.

B)

Add a sinusoid of frequency 220 Hz to the pink noise signal with
power equal to that of the pink noise (i.e. p = std(pink noise signal) and
sig = sqrt(2)*p*sin()) repeat part (A). Show the plots and explain the
differences observed with the plots from part (A). Does this make
sense since the sine function is not random?


Homework (6)

a)

Create a tone at 450 Hz with sampling rate 8000, amplitude
0.707, and duration 3 seconds. Add white noise (use the
randn

function) with a signal
-
to
-
noise ratio of 12 dB. Use
the
fir1

command to design a 30 and 120
th

order band
-
pass
filter from 400 Hz to 506 Hz. Plot the magnitude response
of the filters. Use
filter

to filter the signal and listen to the
sound before and after filtering. Plot the signal spectral
magnitudes before and after filtering (use the
fft

function).
Describe the differences you hear between the signals for
before and after filtering and compare the before and after
filtering spectra.

b)

Repeat part (a) comparing a 5
th

order elliptical filter with
passband ripple of .5 dB and stopband ripple of 30 dB and a
5
th

order Butterworth filter.

Homework(7)


Use Sinc function interpolation to upsample by a factor of
50 the band
-
limited signal is given in terms of its samples
below. Assume signal is bandlimited to 0.5 Hz. The
sampling rate (
B
) is 1 Hz.




F = [0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0]



Note:


a) Matlab has a sinc function (see help sinc)


b) A step edge is not bandlimited so the interpolated
bandlimited function will have overshoot and undershoot
between the given samples (i.e. it will not look like a
sharper step).