Department of Electronics & Telecommunication Engineering, MMCT
Introduction to Matlab,
Signal Processing & Speech
Signal Processing
Outline:
What is Matlab?
Matlab Screen
Variables, array, matrix, indexing
Operators (Arithmetic, relational, logical )
Using of M-File
Writing User Defined Functions
Signal Processing
Speech Signal Processing
What is Matlab?
Matlab is basically a high level language
which has many specialized toolboxes for
making things easier for us
How high?
Matlab
High Level
Languages such as
C, Pascal etc.
Assembly
What are we interested in?
Matlab is too broad for our purpose.
The features we are going to require is
Series of
Matlab
commands
Matlab
m-files
functions
Input
Output
capability
Command
Line
mat-files
Matlab Screen
Command Window
type commands
Current Directory
View folders and m-files
Workspace
View program variables
Double click on a variable
to see it in the Array Editor
Command History
view past commands
save a whole session
using diary
Variables
No need for types. i.e.,
int a;
double b;
float c;
All variables are created with double precision unless
specified and they are matrices.
Example:
>>x=5;
>>x1=2;
After these statements, the variables are 1x1 matrices
with double precision
Array, Matrix
a vector
x = [1 2 5 1]
x =
1
a matrix
x = [1 2 3; 5 1 4; 3 2 -1]
x =
1
5
3
transpose
2
1
2
3
4
-1
y = x
y =
1
2
5
1
Long Array, Matrix
t =1:10
t =
2
3
4
k =2:-0.5:-1
k =
1.5
0.5
-0.5
= [1:4; 5:8]
x =
1
5
2
6
3
7
4
8
-1
10
Matrix Index
The matrix indices begin from 1 (not 0 (as in C))
The matrix indices must be positive integer
Given:
A(-2), A(0)
Error: ??? Subscript indices must either be real positive integers or logicals.
A(4,2)
Error: ??? Index exceeds matrix dimensions.
Concatenation of Matrices
x = [1 2], y = [4 5], z=[ 0 0]
A = [ x y]
1
B = [x ; y]
1 2
4 5
C = [x y ;z]
Error:
??? Error using ==> vertcat CAT arguments dimensions are not consistent.
Operators (arithmetic)
+
*
/
^
addition
subtraction
multiplication
division
power
complex conjugate transpose
Matrices Operations
Given A and B:
Addition
Subtraction
Product
Transpose
Operators (Element by
Element)
.* element-by-element multiplication
./ element-by-element division
.^ element-by-element power
The use of . Element
A =Operation
[1 2 3; 5 1 4; 3 2 -1]
A=
1
5
3
x = A(1,:)
x=
2
1
2
3
4
-1
y = A(3 ,:)
y=
1 2 3
3 4 -1
b = x .* y
c=x./y
d = x .^2
b=
c=
0.33 0.5 -3
d=
3 8 -3
K= x^2
Erorr:
??? Error using ==> mpower Matrix must be square.
B=x*y
Erorr:
??? Error using ==> mtimes Inner matrix dimensions must agree.
General Commands
On-line help
help
helpwin
helpdesk
help topic
lookfor string
demo
Directory Information & Termination
pwd
cd
dir
ls
Path
editpath
copyfile
Workspace information
who
whos
what
clear
clf
clc
home
General Information & Termination
computer
clock
date
ver
^c
quit or exit
Arithmetic Operations
Compute following
5
1.
2
5
(2 - 1)
Result
1.0323
2
2 Area= r
1/3
take r=
3. Sin /6 0.5
x
4.
= 17
for
x 2.57
- 1
0.6781
Basic Task: Plot the function
sin(x) between 0x4
Create an x-array of 100 samples between 0 and
4.
>>x=linspace(0,4*pi,100);
Calculate sin(.) of the x-array
1
0.8
0.6
>>y=sin(x);
Plot the y-array
>>plot(y)
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
10
20
30
40
50
60
70
80
90
100
Plot the function e-x/3sin(x)
between 0x4
Create an x-array of 100 samples between 0
and 4.
>>x=linspace(0,4*pi,100);
Calculate sin(.) of the x-array
>>y=sin(x);
Calculate e-x/3 of the x-array
>>y1=exp(-x/3);
Multiply the arrays y and y1
>>y2=y*y1;
Plot the function e-x/3sin(x)
between 0x4
Multiply the arrays y and y1 correctly
>>y2=y.*y1;
Plot the y2-array
>>plot(y2)
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
-0.1
-0.2
-0.3
10
20
30
40
50
60
70
80
90
100
Display Facilities
0.7
0.6
plot(.)
0.5
0.4
0.3
Example:
>>x=linspace(0,4*pi,100);
>>y=sin(x);
>>plot(y)
>>plot(x,y)
stem(.)
0.2
0.1
0
-0.1
-0.2
-0.3
10
20
30
40
50
60
70
80
90
100
10
20
30
40
50
60
70
80
90
100
0.7
0.6
0.5
0.4
0.3
Example:
>>stem(y)
>>stem(x,y)
0.2
0.1
0
-0.1
-0.2
-0.3
Display Facilities
title(.)
>>title(This is the sinus function)
xlabel(.)
>>xlabel(x (secs))
0.8
ylabel(.)
0.6
0.4
0.2
sin(x)
This is the sinus function
0
-0.2
-0.4
-0.6
>>ylabel(sin(x))
-0.8
-1
10
20
30
40
50
60
x (secs)
70
80
90
100
Operators (relational,
logical)
== Equal to
~= Not equal to
< Strictly smaller
> Strictly greater
<= Smaller than or equal to
>= Greater than equal to
& And operator
| Or operator
Flow Control
if
for
while
break
.
Control Structures
If Statement Syntax
if (Condition_1)
Matlab Commands
elseif (Condition_2)
Matlab Commands
elseif (Condition_3)
Matlab Commands
else
Matlab Commands
end
Some Dummy Examples
if ((a>3) & (b==5))
Some Matlab Commands;
end
if (a<3)
Some Matlab Commands;
elseif (b~=5)
Some Matlab Commands;
end
if (a<3)
Some Matlab Commands;
else
Some Matlab Commands;
end
Control Structures
Some Dummy Examples
For loop syntax
for i=Index_Array
Matlab Commands
end
for i=1:100
Some Matlab Commands;
end
for j=1:3:200
Some Matlab Commands;
end
for m=13:-0.2:-21
Some Matlab Commands;
end
for k=[0.1 0.3 -13 12 7 -9.3]
Some Matlab Commands;
end
Control Structures
While Loop Syntax
while (condition)
Matlab Commands
end
Dummy Example
while ((a>3) & (b==5))
Some Matlab Commands;
end
Use of M-File
Click to create
a new M-File
Extension .m
A text file containing script or function or program to run
Use of M-File
Save file as Denem430.m
If you include ; at the
end of each statement,
result will not be shown
immediately
Writing User Defined Functions
Functions are m-files which can be executed by
specifying some inputs and supply some desired outputs.
The code telling the Matlab that an m-file is actually a
function is
function out1=functionname(in1)
function out1=functionname(in1,in2,in3)
function [out1,out2]=functionname(in1,in2)
You should write this command at the beginning of the mfile and you should save the m-file with a file name same
as the function name
Writing User Defined Functions
Another function which takes an input array and returns the sum and product
of its elements as outputs
The function sumprod(.) can be called from command window or an m-file as
Useful Commands
The two commands used most by Matlab
users are
>>help functionname
>>lookfor keyword
Signal Processing Examples
To generate a cosine signal.
clc;
t=0:.01: pi;
y=cos (2*pi*t);
subplot (2,1,1);
plot (t,y);
subplot(2,1,2);
stem(t,y);
display('the resultant is');
Program Output:
1
0.5
0
-0.5
-1
0.5
1.5
2.5
3.5
0.5
1.5
2.5
3.5
1
0.5
0
-0.5
-1
To generate the exponential
signal
clc;
n= input('enter the length of exponential
signal');
t= 0:n;
a= input('enter the a value');
z= exp(a*t);
subplot(2,2,1);
stem(z,t);
display('the result length of exponential
signal');z
Output:
Enter the length of exponential signal 4
Enter the a value 1
The result length of exponential signal
z=
1.0000 0.3679 0.1353 0.0498 0.0183
Writing User Defined Functions
Examples
Write a function : out=squarer (A, ind)
Which takes the square of the input matrix if the input
indicator is equal to 1
And takes the element by element square of the input
matrix if the input indicator is equal to 2
Same Name
DFT
Discrete Fourier transform (DFT) converts
a finite list of equally spaced samples of a
function into the list of coefficients of a finite
combination of complex sinusoids, ordered
by theirfrequencies, that has those same
sample values. It can be said to convert the
sampled function from its original domain
(often time or position along a line) to
the frequency domain.
DFT
In digital signal processing, the function is
any quantity orsignal that varies over time,
such as the pressure of a sound wave,
a radio signal, or daily temperature readings,
sampled over a finite time interval (often
defined by a window function). In image
processing, the samples can be the values
of pixelsalong a row or column of a raster
image
DFT
The sequence of N complex numbers is
transformed into an N-periodic sequence of
complex numbers according to the DFT
formula:
N-1
-j(2/N)nk
F[n]= f[k] e
(n=0; N-1)
k=0
Computing
DiscreteFourierTransform
clc;
xn=[2,1,3,1];
n=50;
l=length(xn);
x1=[xn,zeros(1,n-l)];
for k=0:1:n-1;
for g=0:1:n-1;
p=exp(-i*2*pi*g*k/n);
x2(k+1,g+1)=p;
end
end
xk=x1*x2;
subplot(2,2,1);
y= abs(xk);
stem(y);
z=angle(xk);
subplot(2,2,2);
stem(z);
disp('the resultant signal is');y
Output
8
-2
20
40
60
-4
20
40
60
Matlab Functionality in
Speech Signal Processing
Basic Functionality
read a speech file (i.e., open a .wav speech file and read the speech sample into a
MATLAB array)
write a speech file (i.e., write a MATLAB array of speech samples into a .wav speech file)
play a MATLAB array of speech samples as an audio file
plot a speech file (MATLAB array) as a waveform using a strips plot format
plot a speech file (MATLAB array) as one or more 4 line plot(s)
convert the sampling rate associated with a speech file (MATLAB array) to a
different sampling rate
high pass filter a speech file (MATLAB array) to eliminate hum and low frequency noise
plot a frame of speech and its associated spectral log magnitude
plot a spectrogram of a speech file (MATLAB array)
Recording of speech signal
clc;
close all;
%please say your names to record it
For i = 1:4
Fs=22050;
file = sprintf('%s%d.wav','g',i);
input('You have 2 seconds to say your
record--> ');
y = wavrecord(2*Fs,Fs);
sound(y,Fs);
wavwrite(y,Fs,file);
end
name. Press enter when ready to
Naming the Speech File
for i = 1:4
strI=int2str(i);
wavName='g';
temp=strcat(wavName,strI);
%get the speech signal sample
[x,fs]=wavread(temp,[1 44100]);
ms1=fs/1000;
ms20=fs/50;
% maximum speech Fx at 1000Hz
% minimum speech Fx at 50Hz
% hamming window
Y=fft(x.*hamming(length(x)));
Waveform Plot
% plot waveform
t=(0:length(x)-1)/fs;
% times of sampling instants
figure;
name1='waveform of '
title1=strcat({name1},{temp});
plot(t,x);
Grid
title(title1);
xlabel('Time (s)');
ylabel('Amplitude');
end
Plotting the spectrum
% plot spectrum of bottom 5000Hz
hz5000=5000*length(Y)/fs;
f=(0:hz5000)*fs/length(Y);
figure;
name2=' spectrum of '
title2=strcat({name2},{temp});
plot(f,20*log10(abs(Y(1:length(f)))+eps));
grid;
title(title2);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
Cepstral Analysis to find fundamental
frequency
% cepstrum is DFT of log spectrum
C=fft(log(abs(Y)+eps));
% plot between 1ms (=1000Hz) and 20ms (=50Hz)
q=(ms1:ms20)/fs;
figure;
name3=' cepstrum of '
title3=strcat({name3},{temp});
plot(q,abs(C(ms1:ms20)));
grid;
title(title3);
xlabel('Quefrency (s)');
ylabel('Amplitude');
%search for the index of the peak in the cepstrum between 1 and 20ms, and then convert back to hertz
[c,fx]=max(abs(C(ms1:ms20)));
fprintf('fundamentalfrequency=%gHz\n',fs/(ms1+fx-1));
%fundamental frequency estimation time domain
% get a section of vowel
[x,fs]=wavread(temp,[1 44100]);
ms20=fs/50;
% minimum speech Fx at 50Hz
Fundamental Frequency (Speech
Feature)
%search for the index of the peak in the cepstrum between 1 and 20ms, and then convert back to
hertz
[c,fx]=max(abs(C(ms1:ms20)));
fprintf('fundamentalfrequency=%gHz\n',fs/(ms1+fx-1));
Thank You
