Lab 2 – Elementary Signals Objective: To study and plot the models of some elementary signals.
Theory: A signal can be defined as a function fu nction of one or more variables that conveys information about some physical phenomenon. Signals are modeled by mathematical functions. In the study of signals and systems, several elementary signals feature prominently. These include the step function, the impulse function, exponential signals, and sinusoidal signals. Not only that these signals occur naturally in a wide range of physical systems, they also serve as building blocks for constructing more complex signals. In this lab, we will study such signals and learn to model these signals using Matlab. We will limit signals to having one independent variable only. For continuous-time signals, the independent variable will be time t, and for discrete-time signals, the independent variable will be the sample number n. We will be modeling the following signals: The Unit Step Function The unit step function in continuous-time is defined as:
0, t 0 . u (t ) 1, t 0 Notice that the function has a discontinuity at t 0 . In discrete-time, the unit step function is defined as:
0, n 0 u[n] . 1 , 0 n The Impulse Function The continuous-time unit-impulse function is defined as: a s:
1, t 0 . 0, 0 t
(t )
The unit impulse signal in discrete-time is given as:
1, n 0 . 0, 0 n
[n]
It should be noted that in case of continuous-time impulse, the value of (t ) specifies the weight rather than the amplitude of the impulse. The Real Exponential A continuous-time real exponential signal is defined as: x(t ) Ae , at
where A is the amplitude of the signal at t 0 and while for 1 , the signal grows exponentially.
is a real number. For 1 , the signal decays
The discrete-time real exponential signal is expressed as: x[n] Ae , n
where A and
have same meanings as in the previous equation. .
The Sinusoidal Signal In its most general form, the continuous-time sinusoidal signal can be given as: x(t ) A cos(0t )
where A is the amplitude, 0 is the angular frequency in radians per second, and is the phase angle in radians. The continuous-time sinusoidal signal is a periodic signal si gnal with time period given by T 2 / 0 . The discrete-time sinusoidal signal can be expressed as: x[n] A cos(n )
If n is taken to be dimensionless, then both and have units of radians. It must be noted that a discrete-time sinusoidal signal may or may not be periodic. For it to be periodic with a period of N samples, must be a rational multiple of 2 .
Modeling Elementary Signals Now, we will model the above signals using Matlab. In general, these signal models are defined for all t and all n , however, we will be modeling the signals within a small time-span around n 0 or t 0 . It should also be noted that Matlab works on discrete-time signals. Hence, we will only be approximating continuous-time signals by keeping a very small sampling interval and interpolating the samples using plot plot function of Matlab. Type and run the following programs in you m-file editor and run the programs.
%% Lab 2.1 - Generating Unit Impulse Function %% Initializing close all clear all clc N=10;% N=10;% Number of sample on both sides of n = 0 n=[-N:1:N]; %% Generating the Impulse imp=zeros(1,length(n));%Generating imp=zeros(1,length(n)); %Generating a vector of zeros imp(n==0)=1;%Setting imp(n==0)=1;%Setting the sample at n=0 to 1 %% Plotting the function stem(n,imp,'LineWidth' stem(n,imp,'LineWidth',2); ,2); axis([n(1) n(length(n)) min(imp)-0.1 max(imp)+0.1]); xlabel('Sample xlabel('Sample Number'); Number'); ylabel('Magnitude' ylabel('Magnitude'); ); title('Discrete-time title('Discrete-time Unit Impulse Function'); Function'); grid on on; ;
The output of the program can be seen in Figure 2.1. Discrete-time Unit Impulse Function
1
0.8
e d u t i n g a M
0.6
0.4
0.2
0
-10
-8
-6
-4
-2
0 2 Sample Number
Figure 2. 1. Output of Lab 2.1
4
6
8
10
%% Lab 2.2 - Generating Unit Step Function %% Initializing close all clear all clc %% Discrete-time N=10;% N=10;% Number of sample on both sides of n = 0 n=[-N:1:N]; %% Generating the Step Function DiscreteStep=zeros(1,length(n));%Generating DiscreteStep=zeros(1,length(n)); %Generating a vector of zeros DiscreteStep(n>=0)=1;%Setting DiscreteStep(n>=0)=1; %Setting the sample at n>=0 to 1 %% Plotting the function figure subplot(211); stem(n,DiscreteStep,'LineWidth' stem(n,DiscreteStep, 'LineWidth',2); ,2); axis([n(1) n(length(n)) min(DiscreteStep)-0.1 max(DiscreteStep)+0.1]); xlabel('Sample xlabel('Sample Number'); Number'); ylabel('Magnitude' ylabel('Magnitude'); ); title('Discrete-time title('Discrete-time Unit Step Function'); Function'); grid on on; ; %% Continuous-time T=2;% T=2;% Time on both sides of t = 0 t=[-T:0.0011:T]; % %% Generating the Step Function ContStep=zeros(1,length(t));%Generating ContStep=zeros(1,length(t)); %Generating a vector of zeros ContStep(t>0)=1;%Setting ContStep(t>0)=1; %Setting the sample at n=0 to 1 %% Plotting the function subplot(212) plot(t,ContStep,'LineWidth' plot(t,ContStep, 'LineWidth',2); ,2); axis([t(1) t(length(t)) min(ContStep)-0.1 max(ContStep)+0.1]); xlabel('Time' xlabel('Time'); ); ylabel('Magnitude' ylabel('Magnitude'); ); title('Continuous-time title('Continuous-time Unit Step Function'); Function'); grid on on; ;
Notice the sampling interval of 0.0011 seconds to generate the time vector. This sampling interval avoids the signal to be sampled at t = 0 at which the function is undefined. The output of the program is shown in Figure 2.2.
Discrete-time Unit Step Function 1 0.8 e d u 0.6 t i n g a 0.4 M
0.2 0 -10
-8
-6
-4
-2
0 2 Sample Number
4
6
8
Continuous-time Unit Step Function 1 0.8 e d u 0.6 t i n g a 0.4 M
0.2 0 -2
-1.5
-1
-0.5
0 Time
0.5
Figure 2. 2. Output of Lab 2.2 %% Lab 2.3 - Generating Exponential Signals %% Initializing close all clear all clc %% Discrete-time N=10;% N=10;% Number of sample on both sides of n = 0 n=[-N:1:N]; alpha = 0.1; A=1; DiscreteExp1=A*exp(alpha*n); DiscreteExp2=A*exp(-alpha*n);
1
1.5
10
%% Plotting the signals figure subplot(221); stem(n,DiscreteExp1,'LineWidth',2); axis([n(1) n(length(n)) min(DiscreteExp1)-0.1 max(DiscreteExp1)+0.1]); xlabel('Sample Number'); ylabel('Magnitude'); grid on subplot(222); stem(n,DiscreteExp2,'LineWidth',2); axis([n(1) n(length(n)) min(DiscreteExp2)-0.1 max(DiscreteExp2)+0.1]); xlabel('Sample Number'); ylabel('Magnitude'); grid on; %% Continuous-time T=2;% Number of sample on both sides of n = 0 t=[-T:0.001:T]; alpha = 0.5; A=1; ContExp1=A*exp(alpha*t); ContExp2=A*exp(-alpha*t); %% Plotting the signals subplot(223); plot(t,ContExp1,'LineWidth',2); axis([t(1) t(length(t)) min(ContExp1)-0.1 max(ContExp1)+0.1]); xlabel('Time'); ylabel('Magnitude'); grid on subplot(224); plot(t,ContExp2,'LineWidth',2); axis([t(1) t(length(t)) min(ContExp2)-0.1 max(ContExp2)+0.1]); xlabel('Time'); ylabel('Magnitude'); grid on;
The output can be seen in Figure 2.3.
2.5
2.5
2 e d u t i n 1.5 g a M
2 e d u t i n 1.5 g a M
0.5
0.5
1
-10
1
-5
0 5 Sample Number
10
-10
2.5
2.5
2 e d u t i n 1.5 g a M
2 e d u t i n 1.5 g a M
0.5
0.5
1
-2
-5
0 5 Sample Number
10
1
-1
0 Time
1
2
-2
Figure 2. 3. Output of Lab 2.3
%% Lab 2.4 - Generating Sinusoidal Signals %% Initializing close all clear all clc %% Continuous-time T=2; t=[-T:0.001:T]; A=2; omega=2*pi; phi=pi/4; x=A*sin(omega*t+phi);
-1
0 Time
1
2
%% Plotting subplot(211); plot(t,x,'LineWidth',2); title('Sinusoidal Signal'); axis([-T T min(x)-0.1 max(x)+0.1]); xlabel('Sample Number'); ylabel('Magnitude'); grid on %% Discrete-time N=20; n=[-N:1:N]; A=2; omega=pi/8; phi=pi/4; y=A*sin(omega*n+phi); %% Plotting subplot(212); stem(n,y,'LineWidth',2); axis([-N N min(y)-0.1 max(y)+0.1]); xlabel('Time'); ylabel('Magnitude') grid on
The output of the program can be seen in Figure 2.4. Sinusoidal Signal 2
e d u t i n g a M
1 0 -1 -2 -2
-1.5
-1
-0.5
-15
-10
-5
0 0.5 Sample Number
1
1.5
2
10
15
20
2
e d u t i n g a M
1 0 -1 -2 -20
0 Time
5
Figure 2. 4. Output of Lab 2.4
Exercises Write Matlab code to generate and plot the following signals: 1. x(t ) Ae t where A 1 . Try different values of and observe the effect of changing the values.
2. For the continuous-time signal in Lab 2.4, what is the time delay corresponding to the phase angle of / 4 radians. Change the value of to 0 and / 2 and comment on the result.
3. For the discrete-time signal in Lab 2.4, change the values of to / 2 and / 3 and comment on the periodicity of the signal.
Lab 3 – Basic Signal Operations Objective: To perform basic signal operations on the elementary signals.
Theory: Basic operations on signals can be classified into two classes: 1. Operations performed on dependent variables: a. Amplitude scaling b. Addition c. Multiplication 2. Operations performed on independent variables: a. Time scaling b. Reflection c. Time shifting Amplitude Scaling For a continuous-time signal x(t ) , the signal y (t ) resulting from amplitude scaling applied to x(t ) is given by: y(t ) cx(t )
where c is the scaling factor. Likewise, for a discrete-time signal x[n ] , the signal y[n] resulting from amplitude scaling applied to x[n] is given by: y[n] cx[n]
Signal Addition
Let x1 (t ) and x2 (t ) denote a pair of discrete-time signals, then the signal y(t ) obtained by the addition of x1 (t ) and x2 (t ) is defined by: y (t ) x1 (t ) x2 (t )
Likewise, addition of two discrete-time signals is given by: y[n] x1[n] x2 [n]
Signal Multiplication
Let x1 (t ) and x2 (t ) denote a pair of discrete-time signals, then the signal y(t ) obtained by the multiplication of x1 (t ) and x2 (t ) is defined by: y(t ) x1 (t )x2 (t )
Similarly, in discrete-time, multiplication of two signals can be expressed as: y[n] x1[n]x2 [n ]
This means that for each n ¸ the value of y[ n] is given by the product of the corresponding values of x1[ n] and x2 [ n ] .
Reflection Given a signal x(t ) , its time reversed or reflected version is given by y(t) x(t ) . The signal y(t )
represents a reflected version of x(t ) about
t 0
.
Time Scaling Let x(t ) denotes a continuous-time signal. Then the signal y(t ) obtained by scaling the independent
variable, time t , by a factor a is defined as y(t ) x(at )
If a > 1, y(t ) is a compressed version of x(t ) . If 0 < a < 1, then y(t ) is an extended version of x(t ) . Time Shifting Given a signal x(t ) , its time-shifted version is defined by: y(t ) x(t t 0 )
where t 0 is the time-shift. For t 0 0 , the waveform of y(t ) is obtained by shifting x(t ) towards the right, relative to the time-axis. If t 0 0 , x(t ) is shifted to the left.
Basic Signal Operations using Matlab Run the following Matlab programs. The output of Lab 3.1 can be seen in Figure 3.1, while the output of Lab 3.2 is shown in Figure 3.2.
%% Lab 3.1 - Amplitude scaling %% Initializing close all clear all clc %% Generating Exponential Signal N=10;% Number of sample on both sides of n = 0 n=[-N:1:N]; alpha = 0.1; A=1; DiscreteExp=A*exp(alpha*n); %% Amplitude Scaling ScaledExp1=2*DiscreteExp; ScaledExp2=0.3*DiscreteExp; %% Plotting the signals figure subplot(131); stem(n,DiscreteExp,'LineWidth',2); axis([n(1) n(length(n)) min(DiscreteExp)-0.1 max(DiscreteExp)+0.1]); xlabel('Sample Number'); ylabel('Magnitude'); grid on subplot(132); stem(n,ScaledExp1,'LineWidth',2); axis([n(1) n(length(n)) min(ScaledExp1)-0.1 max(ScaledExp1)+0.1]); xlabel('Sample Number'); ylabel('Magnitude'); grid on; subplot(133); stem(n,ScaledExp2,'LineWidth',2); axis([n(1) n(length(n)) min(ScaledExp2)-0.1 max(ScaledExp2)+0.1]); xlabel('Sample Number'); ylabel('Magnitude'); grid on;
%% Lab 3.2 - Signal Addition and Multiplication %% Initializing close all clear all clc N=10;% Number of sample on both sides of n = 0 n=[-N:1:N]; %% Generating Exponential Signal alpha = 0.1; A=1; DiscreteExp=A*exp(alpha*n); %% Generating Sinusoidal Signal A=2; omega=pi/4; phi=0; SineWave=A*sin(omega*n+phi); %% Adding two signals y=DiscreteExp+SineWave; %% Multiplying two signals %multiplication %% Plotting the signals figure subplot(221); stem(n,DiscreteExp,'LineWidth',2); axis([n(1) n(length(n)) min(DiscreteExp)-0.1 max(DiscreteExp)+0.1]); xlabel('Sample Number'); ylabel('Magnitude'); grid on title('Signal 1'); subplot(222); stem(n,SineWave,'LineWidth',2); axis([n(1) n(length(n)) min(SineWave)-0.1 max(SineWave)+0.1]); xlabel('Sample Number'); ylabel('Magnitude'); grid on; title('Signal 2'); subplot(223); stem(n,y,'LineWidth',2); axis([n(1) n(length(n)) min(y)-0.1 max(y)+0.1]); xlabel('Sample Number'); ylabel('Magnitude'); grid on; title('Signal1 + Signal2')
subplot(224); stem(n,z,'LineWidth',2); axis([n(1) n(length(n)) min(z)-0.1 max(z)+0.1]); xlabel('Sample Number'); ylabel('Magnitude'); grid on; title('Signal1 x Signal2'); 5.5
0.9
5
2.5
0.8
4.5
4
2
e d u t i n g a M
0.7
0.6 e d u t i n g a M
1.5
3.5
e d u t i n g a M
3
0.5
0.4
2.5 0.3 1
2 0.2 1.5
0.5
-10
0.1
1 0 Sample Number
10
-10
0 Sample Number
Figure 3. 1. Output for Lab 3.1
10
-10
0 Sample Number
10
Lab 4 – Convolution Objective: To perform convolution of continuous-time signals.
Theory: Is a mathematical way of combining two signals to achieve a third, modified signal. The signal we record seems to respond well to being treated as a series of signals superimposed upon each other that is seismic signals seem to respond convolutionally. The process of DECONVOLUTION is the reversal of the convolution process.
Mathematically, a convolution is defined as the integral over all space of one function at x times another function at u-x. The integration is taken over the variable x (which may be a 1D or 3D variable), typically from minus infinity to infinity over all the dimensions. So the convolution is a function of a new variable u, as shown in the following equations. The cross in a circle is used to indicate the convolution operation.
Note that it doesn't matter which function you take first, i.e. the convolution operation is commutative.
Basic Example Let us look at a basic continuous-time convolution example to help express some of the important ideas. We will convolve together two square pulses, x(t ) and
h(t ), as shown in Figure 1
(a)
(b)
Figure 1: Two basic signals that we will convolve together.
Reflect and Shift Now we will take one of the functions and reflect it around the y-axis. Then we must shift the function, such that the origin, the point of the function that was originally on the origin, is labeled as point t . This step is shown in Figure 2,
h(t −τ ).
(a) Reflected square pulse.
(b) Reflected and shifted square pulse.
Figure 2: h(−τ ) and h(t −τ ).
Note that in Figure 2 τ is the 1st axis variable while
t is a constant (in this figure). Since convolution is commutative
it will never matter which function is reflected and shifted; however, as the functions become more complicated reflecting and shifting the "right one" will often make the problem much easier.
Regions of Integration
∫∞−∞ x(τ )h(t −τ )dτ . The value of the function y at time t is given by the amount of overlap(to be precise the integral of the overlapping region) between h(t −τ ) and x(τ ). We start out with the convolution integral, y(t )=
Next, we want to look at the functions and divide the span of the functions into different limits of integration. These
−
different regions can be understood by thinking about how we slide h(t τ ) over x(τ ), see Figure 3.
(a) No overlap.
(b) h(t −τ ) on its way "into" x(τ )
(c) h(t −τ ) on its way "out of" x(τ )
(d) No overlap.
Figure 3: Figures to help understand the regions of intergration
In this case we will have the following four regions. Compare these limits of integration to the four illustrations
−
of h(t τ ) and x(τ ) in Figure 3. Four Limits of Integration 1. 2. 3.
t <0 0≤t <1 1≤t <2
4.
t ≥2
Using the Convolution Integral
−
Finally we are ready for a little math. Using the convolution integral, let us integrate the product of x(τ )h(t τ ). For our first and fourth region this will be trivial as it will always be 0. The second region, 0 following math:
≤t <1, will require the
y(t )==∫t 01dτ=t
(1)
≤t <2, is solved in much the same manner. Take note of the changes in our integration though. As we move h(t −τ ) across our other function, the left-hand edge of the function, t −1, becomes our lowlimit for The third region, 1
the integral. This is shown through our convolution integral as
y(t )===∫1t −11dτ=1−(t −1)=2−t (2)
The above formulas show the method for calculating convolution; however, do not let the simplicity of this example confuse you when you work on other problems. The method will be the same, you will just have to deal with more math in more complicated integrals. Note that the value of y(t ) at all time is given by the integral of the overlapping functions. In this example y for a given t equals the gray area in the plots in Figure 3.
Convolution Results (3)
Now that we have found the resulting function for each of the four regions, we can combine them together and graph the convolution of x(t )∗h(t ).
Figure 4: Shows the system's output in response to
the input, x(t ).
Matlab code % %
Description:
This m-file convolves two signals and plots the result.
x = [0.5 0.5 0.5]; h = [3.0 2.0 1.0];
% define input signal x[n] % define unit-pulse response h[n]
y = conv(x,h);
% compute output y[n] via convolution
n = 0:(length(y)-1);
% create n-vector for plotting y[n]
stem(n,y); % plot y[n] grid; xlabel('n'); ylabel('y[n]'); title('Output of System via Convolution');
Output of System via Convolution 3
2.5
2
] n [ y
1.5
1
0.5
0
0
0.5
1
1.5
2 n
2.5
3
3.5
4
clear all x=[2 3 4 1 2 3 7 8]; t=1:10; h=exp(-t); sec y=conv(x,h); plot(y) grid on xlabel('time') ylabel('output')
%input signal to the system %System response of the system defined from 1 to 10
Convolution 4.5 4 3.5 3 t u p t u o
2.5 2 1.5 1 0.5 0
0
2
4
6
8
10
12
14
16
time
Exercise; αt
Q1) Convolve the A е
and the unit step function when е has the decreasing curve.
Q2) Prove the output of the convolution integral example in the matlab.
18
Lab 5 – Constructing periodic signals from harmonically related sinusoidal signals
Objectives To generate periodic signals by combining harmonically related sinusoidal signals
Theory A signal is periodic if, for some positive value of T, x(t T ) x(t ) for all t. The fundamental period of x(t ) is the minimum positive nonzero value of t for which the above equation is satisfied.
The sinusoidal signal x(t ) A cos(0t ) and the complex exponential signal Ae j ( t ) are both 0
periodic with fundamental frequency 0 and fundamental period T 2 / 0 . Associated with the complex exponential signal is the set of harmonically related complex exponentials given as: 1 2,... k (t ) e jk t e jk (2 / T )t , k 0, 0
A linear combination of harmonically related complex exponentials of the form
x(t )
ak e
k
jk0t
ak e jk (2 / T )t 0
k
is also periodic with period T. In the above equation, the term for k = 0 is a constant, while the terms for both k = 1 and k = -1have fundamental frequency equal to 0 and are referred to as fundamental components or the first harmonic components. In general, the components for k = N and k = -N are called Nth harmonic components. x(t )
Let
3
a e
jk 2 t
k
k 3
where
a0 1 , a1
a1 1/ 4 , a2 a2 1/ 2 , a3 a3 1/ 3 .
The above equation can be re-written by collecting the harmonic components with same fundamental frequencies as follows:
x(t ) 1
1
e 4
j 2 t
e j 2 t
1
e 2
j 4 t
e j 4 t
1
e 3
j6 t
e j6 t
Using Euler’s theorem: x(t ) 1
1 2
cos 2 t cos 4 t
2 3
cos6 t
In the following code, we will add these harmonically related components one by one to see the effect of adding each component on the shape of the resulting signal.
MATLAB CODE for constructing periodic signals from harmonically related signals %% Lab 5.1 %% Initializing close all clear all clc t=[-5:0.01:5]; % time indexing vector h0=1 h1=0.5*cos(2*pi*t); % returns the cosine value for all values of t h2=cos(4*pi*t); h3=2/3*cos(6*pi*t); figure; subplot(2,2,1); plot(t,h0); xlabel('Time'); ylabel('h0'); subplot(2,2,2); plot(t,h0+h1); xlabel('Time'); ylabel('h0+h1'); subplot(2,2,3); plot(t,h0+h1+h2); xlabel('Time'); ylabel('h0+h1+h2'); subplot(2,2,4); plot(t,h0+h1+h2+h3); xlabel('Time'); ylabel('h0+h1+h2+h3');
Figure 1 Output for Lab 5.1
Type the following code and execute the program:
%% Lab 5.2 %% Initializing close all clear all clc t=[-5:0.01:5]; % time indexing vector x=zeros(1,length(t)); % plotting zeros from 1 till length of t NumOfHarmonics=10; %% L # 5.1 executing a statement a predetermined number of times
for n=1:2:NumOfHarmonics % indicating the starting and ending of loop x=x+4/pi*1/n*sin(n*pi*t/2); end plot(t,x);
Figure 2 Output for lab 5.2
Try different number of harmonics to see the effect of including more harmonics to the signal.
Exercises
Write MATLAB code for generating a signal with the information given below. Vary the number of harmonics being added and comment on the results. 1.
f ( x)
1
2
Assume L = 1.
1
1
n1 n
n t L
sin
t=[-5:0.01:5]; % time indexing vector x=zeros(1,length(t)); % plotting zeros from 1 till length of t NumOfHarmonics=10; %% L # 5.1 executing a statement a predetermined number of times for n=1:1:NumOfHarmonics % indicating the starting and ending of loop x=1/2-(x+1/pi*1/n*sin(n*pi*t)); end figure(1); hold on; grid on; plot(t,x);
______________________________________________________________________________ ______________________________________________________________________________ ______________________________________________________________________________ ______________________________________________________________________________ _______ 2.
f ( x)
8 2
(1)
n1,3,5,...
( n 1) / 2
n
2
n t L .
sin
t=[-5:0.01:5]; % time indexing vector x=zeros(1,length(t)); % plotting zeros from 1 till length of t NumOfHarmonics=10; %% L # 5.1 executing a statement a predetermined number of times for n=1:2:NumOfHarmonics % indicating the starting and ending of loop y= (n-1)/2; y= (-1)^y; x=(x+8/pi^2*(y/n^2)*sin(n*pi*t)); end figure; hold on; grid on; plot(t,x);
3. x(t )
A T
n 1,3,5
2A
n
sin(
n t T
) cos(
2n t ) T
Use A 1, 1, and T 2 . Change the value of to 0.5 and 1.5 and comment on the result.
t=[-5:0.01:5]; % time indexing vector x=zeros(1,length(t)); % plotting zeros from 1 till length of t A=1; T=2; L=1.5; NumOfHarmonics=10; %% L # 5.1 executing a statement a predetermined number of times for n=1:2:NumOfHarmonics % indicating the starting and ending of loop x=.5+(x+2/(n*pi)*sin(n*pi*L*t/2).*cos(n*pi*L*t)); end figure(1); hold on; grid on; plot(t,x);
___________________________________________________________________________ ___________________________________________________________________________ ___________________________________________________________________________ _________________
Lab 6 - Properties of continuous-time Fourier series-I Objective To study and verify various properties of Fourier series f or continuous time signals
Theory In Lab, we studied Fourier series and learnt how periodic signals can be constructed by adding harmonically related complex exponentials. In this and the next lab, we will study various properties of Fourier series for continuous-time periodic signals.
Let x(t ) and y(t ) denote two periodic signals with period T having Fourier series coefficients denoted by ak and bk respectively, i.e. FS
FS
x(t ) ak and y(t ) bk
Then:
Linearity FS
z(t ) Ax(t ) By(t ) ck
Aak Bbk
Time Shifting FS
x(t t0 ) ak e
jk 0t 0
Time Reversal FS
x(t ) a k
Time Scaling
x( t )
x(t )e k
Verifying the properties using MATLAB Type the following code and execute:
j 0t
%% Lab# 6.1 Title: Properties of Continuous Time Fourier Series, Linearity % Generation of 1st signal xt
close all clear all clc t=-1.5:0.005:1.5; % time indexing vector xcos=cos(2*pi*t); % returns cosine for all values of t xt=xcos>0; subplot(2,1,1); plot(t,xt); xlabel('t'); ylabel('x(t)') title('Peridoic Square Wave (T=1, T1=0.250)') set(gca,'ylim',[-0.1 1.1]);% sets the named properties to specified values on the object identified by gca grid on % Generation of 2nd signal yt
T=1; T1=0.125; lenT=T/0.005; ytemp=zeros(1,lenT); % plotting zeros from 1 till value of lenT lenT1=T1/0.005; ytemp(round(lenT/2)-lenT1:round(lenT/2)+lenT1-1)=ones(1,2*lenT1); yt=[ytemp ytemp ytemp 0]; % The last 0 added to make the size of yt equal to length(t) subplot(2,1,2);plot(t,yt);xlabel('t');ylabel('y(t)') title('Periodic Square Wave (T=1, T1=0.125)') set(gca,'ylim',[-0.1 1.1]); % sets the named properties to specified values on the object identified by gca grid z1t=xt+yt; figure; plot(t,z1t); set(gca,'ylim',[-0.1 2.1]); % sets the named properties to specified values on the object identified by gca grid % FS coefficients of periodic square waves
k=-50:50; T1=0.25; ak=sin(k*2*pi*(T1/T))./(k*pi); % returns sine value of angle % Manual correction for a0 -> ak(51) ak(k==0)=2*T1/T; T1=0.125; bk=sin(k*2*pi*(T1/T))./(k*pi); % Manual correction for b0 -> bk(51)
bk(k==0)=2*T1/T; % Application of linearity property of FS ck=ak+bk; % Reconstruction with M=50 w0=2*pi/T; zt=zeros(1,length(t)); for k=-50:50 zt=zt + ck(k+51)*exp(j*k*w0*t); end figure; plot(t,real(zt));grid;xlabel('t');ylabel('z(t)=x(t)+y(t)') title('Reconstruction from ak+bk''s with 101 terms')
%% L # 6.2 Title: Properties of Continuous Time Fourier Series %%Time Shifting %Generation of periodic square wave
close all clear all clc t=-1.5:0.005:1.5; xcos=cos(2*pi*t); xt=xcos>0; % FS coefficients of periodic square wave
k=-50:50; T=1;T1=0.25; ak=sin(k*2*pi*(T1/T))./(k*pi); % Manual correction for a0 -> ak(51) ak(k==0)=2*T1/T; %Amount of time shift t0=0.25; % FS coefficients of the time shifted signal
w0=2*pi/T; bk=ak.*exp(-j*k*w0*t0); % Reconstruction from bk's with 101 terms (M=50)
yt=zeros(1,length(t)); for k=-50:50 yt=yt + bk(k+51)*exp(j*k*w0*t); end figure subplot(2,1,1); plot(t,xt); xlabel('t'); ylabel('x(t)') title('Periodic Square Wave (T=1, T1=0.25)') set(gca,'ylim',[-0.2 1.2]); grid subplot(2,1,2); plot(t,real(yt)) xlabel('t');ylabel('y(t)=x(t-0.5)') title('Reconstruction from bk''s with 101 terms') set(gca,'ylim',[-0.2 1.2]); % sets the named properties to specified values on the object identified by gca grid
Exercise Write MATLAB codes to prove the time reversal and time scaling properties of MATLAB. __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ __________________________________________________________________ ________________________________________
Lab 7 - Properties of continuous-time Fourier series-II Objective To study and verify various properties of Fourier seri es for continuous time signals
Theory This lab is a continuation of Lab 6, where we studied some properties of Fourier series for continuous time signals. In this lab we are going to study to two more properties namely, multiplication and conjugation.
Let x(t ) and y(t ) denote two periodic signals with period T having Fourier series coefficients denoted by ak and bk , respectively. Then Multiplication
FS
x(t ) y(t ) hk
ab
l
k l
l
It can be observed that the right side of the above equation is actually the convolution sum of the Fourier coefficients ak and bk . Conjugation and conjugate symmetry FS
x (t ) a *
MATLAB CODE for verifying the properties Type the following programs and execute:
*
k
%% L # 7.1 Title: Properties of Continuous Time Fourier Series % Multiplication % Generation of 1Hz cosine
close all clear all clc t=-1.5:0.005:1.5; % time indexing vector xt=cos(2*pi*t); % Generation of periodic square wave yt=xt>0; % FS coefficients of periodic square wave k=-50:50; T=1; T1=0.25; ak=sin(k*2*pi*(T1/T))./(k*pi); % Manual correction for a0 -> ak(51) ak(51)=2*T1/T; % FS coefficients of 1Hz cosine (over k=-1..1) bk=zeros(1,3); bk(1)=0.5;bk(3)=0.5; ck=conv(ak,bk); ck(103)=[];ck(1)=[]; % Reconstruction from ck's with 101 terms (M=50)
w0=2*pi/T; zt=zeros(1,length(t)); for k=-50:50 zt=zt + ck(k+51)*exp(j*k*w0*t); end figure(1); set(gcf,'defaultaxesfontsize',8) subplot(3,1,1); plot(t,xt); ylabel('x(t)') title('1 Hz Cosine') set(gca,'ylim',[-1.2 1.2]); grid subplot(3,1,2); plot(t,yt); ylabel('y(t)') title('Periodic Square Wave (T=1, T1=0.25)') set(gca,'ylim',[-1.2 1.2]); grid subplot(3,1,3); plot(t,real(zt)) xlabel('t'); ylabel('z(t)=x(t)*y(t)') title('Reconstruction from ck''s') set(gca,'ylim',[-1.2 1.2]); grid
%% Lab# 7.2 Title: Properties of Continuous Time Fourier Series % Conjugation and Conjugate Symmetry % Generation of 1Hz cosine (Real part of our signal)
close all clear all clc t=-1.5:0.005:1.5; xt=cos(2*pi*t); % Generation of periodic square wave (Imaginary part of our signal)
yt=xt>0; % Our perodic complex valued signal
zt=xt+j*yt; % FS coefficients of 1Hz cosine, xt (over k=-50:50) ak=zeros(1,101); ak(50)=0.5;ak(52)=0.5; % FS coefficients of periodic square wave, yt k=-50:50; T=1;T1=0.25; bk=sin(k*2*pi*(T1/T))./(k*pi); % Manual correction for a0 -> ak(51) bk(51)=2*T1/T; % FS coefficients of zt (using linearty property of FS)
ck=ak+j*bk; % Flip ck's and conjugate dk=conj(fliplr(ck)); % Reconstruction from dk's with 101 terms (M=50) w0=2*pi/T;zrt=zeros(1,length(t)); for k=-50:50 zrt=zrt + dk(k+51)*exp(j*k*w0*t); end figure(1); set(gcf,'defaultaxesfontsize',8) subplot(2,2,1); plot(t,xt); ylabel('x(t)') title('Real part of z(t), 1 Hz Cosine') set(gca,'ylim',[-1.2 1.2]); grid subplot(2,2,2); plot(t,yt); ylabel('y(t)') title('Imag. part of z(t), Periodic Square Wave') set(gca,'ylim',[-1.2 1.2]);grid subplot(2,2,3);plot(t,real(zrt)) xlabel('t'); ylabel('Re[zr(t)]') title('Real part of Reconstruction from dk''s') set(gca,'ylim',[-1.2 1.2]);grid subplot(2,2,4); plot(t,imag(zrt)) xlabel('t'); ylabel('Im[zr(t)]') title('Imag. part of Reconstruction from dk''s') set(gca,'ylim',[-1.2 1.2]);grid
Exercise Write MATLAB programs to prove that:
If x(t ) is real and even, then ak a k ______________________________________________________________________________ ______________________________________________________________________________ ______________________________________________________________________________ ______________________________________________________________________________ ____________________________________________________________________ If x(t ) is real and odd, then ak a k ______________________________________________________________________________ ______________________________________________________________________________ ____________________________________
Lab 8 – introduction to Symbolic Toolbox
Objective: To study the use of Symbolic Toolbox for calculus computation and Fourier series.
Theory: Symbolic computation or algebraic computation or computer algebra relates to algorithms and software for manipulating mathematical expressions and equations in symbolic form, as opposed to manipulating the approximations of specific numerical quantities represented by those symbols. Software applications that perform symbolic calculations are called computer algebra systems .
These systems might be used for symbolic integration or differentiation, substitution of one expression into another, simplification of an expression, etc., for most operations of calculus and, more generally, for every computation with mathematical objects for which algorithms are known. Computer algebra softwares are widely used in many scientific and engineering domains. Symbolic computation is also sometimes referred to as symbolic manipulation, symbolic processing, symbolic mathematics , or symbolic algebra, but these terms also refer to noncomputational manipulation.
Computer Algebra Calculation of Fourier Coefficients A computer algebra system can greatly ease the burden of calculation of the Fourier coefficients of a given function f ( ). t In the case of a function defined "piecewise," we must take care to "split" the integral according to the different intervals of definition of the function. In the paragraphs that follow we illustrate the use of MATLAB (Symbolic toolbox) in deriving the Fourier series
of the period 2square wave function defined on (, ) by
In this case the function is defined by different formulas on two different intervals, so each Fourier coefficient integral from – to must be calculated as the sum of two integrals:
To practice the symbolic derivation of Fourier series in this manner, you can begin by verifying the Fourier series calculated manually. The period 2triangular wave and trapezoidal wave functions illustrated in the figures below have especially interesting Fourier series that we invite you to discover for yourself.
MATLAB code %% Lab 8.1 %% Initializing close all clear all clc syms n t pi %define the cosine coefficients an = (1/pi)*(int(-cos(n*t),-pi,0)+int(cos(n*t),0,pi)) %define the sine coefficients bn = (1/pi)*(int(-sin(n*t),-pi,0)+int(sin(n*t),0,pi)); pretty(bn) %MATLAB does not yet know that n is an integer bn = subs(bn,'(-1)^n','cos(pi*n)'); pretty(bn) %set up a typical Fourier sum FourierSum = (4/pi)*sin(t); for k = 3:2:50 FourierSum = FourierSum+subs((4/pi)*sin(n*t)/n,k,n); end FourierSum %plot its graph ezplot(FourierSum, 3.1416*[-2 4]) grid
Figure 3 (output of Lab 8.1)
Lab 9 – Fourier Transform Objectives
To evaluate the Fourier Transform of signal and plot its real and imaginary part by using MATLAB.
Calculate the inverse Fourier Transform.
Calculate the Power Spectral Density of the define signal.
Theory A signal is periodic if, for some positive value of T, x(t T ) x(t ) for all t. The fundamental period of x(t ) is the minimum positive nonzero value of t for which the above equation is satisfied. The sinusoidal signal x(t ) A cos(0t ) and the complex exponential signal Ae j ( t ) are both periodic with fundamental frequency 0 and fundamental period T 2 / 0 . Associated with the 0
complex exponential signal is the set of harmonically related complex exponentials given as:
k (t ) e jk
0t
e jk (2
, k 0, 1 2,...
/ T )t
A linear combination of harmonically related complex exponentials of the for m
ae
x(t )
jk0 t
k
ak e jk (2 / T )t
k
0
k
is also periodic with period T. In the above equation, the term for k = 0 is a constant, while the terms for both k = 1 and k = -1have fundamental frequency equal to 0 and are referred to as fundamental components or the first harmonic components. In general, the components for k = N and k = -N are called Nth harmonic components. x(t )
Let
3
a e
jk 2 t
k
k 3
Where a0 1 ,
a1
a1
1/ 4 ,
a2
a2
1/ 2 , a3
a3
1/ 3 .
The above equation can be re-written by collecting the harmonic components with same fundamental frequencies as follows: x(t) 1
1
e 4
j 2 t
e j 2 t
1
e 2
j 4 t
e j 4 t
1
e 3
Using Euler’s theorem: x(t ) 1
1 2
cos2 t cos 4 t
2 3
cos6 t
j 6 t
e j6 t
Fourier analysis Fourier analysis follows from Fourier’s theorem, which states that every function can be completely expressed as a sum of sines and cosines of various amplitudes and frequencies. This is a pretty impressive assertion – no matter what the shape of a function, and how little it looks like a sine wave, it can be rewritten as a sum of sines and cosines. The Fourier series tells yo u the amplitude and frequency of the sines and cosines that you should add up to recreate your original function.
Power spectral density The above definitions of energy spectral density require that the Fourier transforms of the signals exist, that is, that the signals are integrable/summable or square-integrable/square-summable. (Note: The integral definition of the Fourier transform is only well-defined when the function is integrable. It is not sufficient for a function to be simply square-integrable. In this case one would need to use the Plancherel theorem.) An often more useful alternative is the power spectral density (PSD), which describes how the power of a signal or time series is distributed with frequency. Here power can be the actual physical power, or more often, for convenience with abstract signals, can be defined as the squared value of the signal, that is, as the actual power dissipated in a load if the signal were a voltage applied across it. This instantaneous power (the mean or expected value of which is the average power) is then given by P(t) = s(t) 2
for a signal s(t). Since a signal with nonzero average power is not square integrable, the Fourier transforms do not exist in this case. Fortunately, the Wiener – Khinchin theorem provides a simple alternative. The PSD is the Fourier transform of the autocorrelation function, R(τ), of the signal if the signal is treated as a wide-sense stationary random process These results are expressed in the mathematical formula,
Matlab Examples
%% Lab 8.1 %% Initializing close all clear all clc N = 8; %% number of sample t = [0:N-1]/N; %% define time f = sin(2*pi*t); %%define function q=fft(f) p = abs(q)/(N/2); %% absolute value of the fft z=ifft(q) subplot(3,2,1) plot(t,f) xlabel('Time on x axis') ylabel('Amplitude') title('Orginal signal') grid on subplot(3,2,2) plot(t,q) xlabel('Time on x axis') ylabel('Fourier Transform ') title('Transformed signal') grid on subplot(3,2,3) plot(t,real(q)) xlabel('Time on x axis') ylabel('Real Part') title('Real Part of the Fourier Transform') grid on subplot(3,2,4) plot(t,p) xlabel('Time on x axis') ylabel('Imiginary Part') title('Imiginary Part of the Fourier Transform') grid on subplot(3,2,5) plot(t,z) xlabel('Time on x axis') ylabel('Ampltiude') title('Original Signal') grid on
Orginal signal 1 e d u t i l p m A
10
0
-0.5
10 t r a P l a e R
-16
Transformed signal
m r o f s 5 n a r T r e 0 i r u o F
0.5
-1
x 10
0 x 10
0.1
0.2
-16
0.3
0.4 0.5 Time on x axis
0.6
0.7
0.8
0.9
-5
0
0.1
0.2
Real Part of the Fourier Transform
0.6
0.7
0.8
0.9
0.7
0.8
0.9
Imiginary Part of the Fourier Transform t r a P y r a 0.5 n i g i m I
0
0.1
0.2
0.3
0.4 0.5 Time on x axis
0.6
0.7
0.8
0.9
0.6
0.7
0.8
0.9
0
0
0.1
0.2
Original Signal 1 0.5 0
-0.5 -1
0.4 0.5 Time on x axis
1
5
-5 0
e d u i t l p m A
0.3
0
0.1
0.2
0.3
0.4 0.5 Time on x axis
Figure 4 Output for Lab 8.1
0.3
0.4 0.5 Time on x axis
0.6
%% Lab 8.2 %% Initializing close all clear all clc N = 10000; %% number of points T = 3.4; %% define time of interval, 3.4 seconds t = [0:N-1]/N; %% define time t = t*T; %% define time in seconds f = sin(2*pi*10*t); %%define function, 10 Hz sine wave p = abs(fft(f))/(N/2); %% absolute value of the fft p = p(1:N/2).^2 %% take the power of positve freq. half freq = [0:N/2-1]/T; %% find the corresponding frequency in Hz semilogy(freq,p); %% plot on semilog scale axis([0 20 0 1]); %% zoom in grid on
0
10
-5
10
-10
10
-15
10
-20
10
-25
10
-30
10
-35
10
0
2
4
6
8
10
12
14
16
18
20
Figure 5 Power spectrum of a 10 Hz sine wave sampled at 10,000 points over 3.4 seconds.
Exercise Q1
Create a MATLAB function that takes a data set and returns the power spectrum. The input arguments to the function should be the data and the actual time spanned by the data. The function should return two vectors, the frequency and the power.
Lab 10 – Applications of Signal & Systems Objective: To design and understand the working of the different parameter s of the filter.
Theory:
Designing the Filter In electronics, computer science and mathematics, a digital filter is a system that performs mathematical operations on a sampled, discrete-time signal to reduce or enhance certain aspects of that signal. This is in contrast to the other major type of electronic filter, the analog filter, which is an electronic circuit operating on continuous- time analog signals. An analog signal may be processed by a digital filter by first being digitized and represented as a sequence of numbers, then manipulated mathematically, and then reconstructed as a new analog signal (see digital signal processing). In an analog filter, the input signal is "directly" manipulated by the circuit. A digital filter system usually consists of an analog-to-digital converter (to sample the input signal), a microprocessor (often a specialized digital signal processor), and a digital-to-analog converter. Software running on the microprocessor can implement the digital filter by performing the necessary mathematical operations on the numbers received from the ADC. In some high performance applications, an FPGA or ASIC is used instead of a general purpose microprocessor. Practice Example
%% Lab 10.1 - design the nth order of the filter %% Initializing close all clear all clc f = [0 .2 .22 .5 .52 1]; %Specify f vector m = [1 1 .4 .4 0 0]; %Specify Amplitude respose n=10; %Specify Order [b a]=yulewalk(n,f,m); %Compute Transfer Function [h w]=freqz(b,a,512); %Determine Filter Response subplot(2,1,1) plot(f/2,m,w/(2*pi),abs(h)) %plot amplitude response grid on xlabel('Normalized frequency r') ylabel('Amplitude') phase=angle(h) %Calculate Phase phase=unwrap(phase) %Unwrap Phase phasedeg=phase*180/pi; %Express Phase into Degree subplot(2,1,2) plot(w/(2*pi),phasedeg)
Sampling and aliasing This is the conversion of a continuous signal into a discrete time signal obtained by taking “samples” of the continuous time signal at discrete time instant. Thus, if xa(t) is the input to the sampler, the output is, xa(nT)=x(n) where t is called sampling interval. Nyquist sampling theorem defines the minimum sampling rate to avoid aliasing. That is minimum sampling rate must be twice that of the highest frequency component of the signal. Aliasing the phenomenon that results in a loss of information when a signal is reconstructed from its sampled signal. In principle, the analog signal can be reconstructed from the samples, provided that the sampling rate is sufficiently high to avoid the problem called “aliasing”. Practice Examples
% Program # 10.2--the sampling of a cosine wave with different sampling % period.
% Illustration of the Sampling Process % in the Time-Domain %% Initializing close all clear all clc clf; t = 0:0.0005:1; f = 13; xa = cos(2*pi*f*t); subplot(2,1,1) plot(t,xa);grid xlabel('Time, msec');ylabel('Amplitude'); title('Continuous-time signal x_{a}(t)'); axis([0 1 -1.2 1.2]) subplot(2,1,2); T = 0.1; % Sampling period n = 0:T:1; xs = cos(2*pi*f*n); k = 0:length(n)-1; stem(k,xs);grid; xlabel('Time index n');ylabel('Amplitude'); title('Discrete-time signal x[n]'); axis([0 (length(n)-1) -1.2 1.2])
% Program P # 10.3--the sampling of a cosine wave with different sampling % period and then reconstruct it highlighting the loss of transmittion.
% Illustration of Aliasing Effect in the Time-Domain %% Initializing close all clear all clc clf; t = 0:0.0005:1; f = 13; xa = cos(2*pi*f*t); subplot(2,1,1) plot(t,xa);grid xlabel('Time, msec');ylabel('Amplitude'); title('Continuous-time signal x_{a}(t)'); axis([0 1 -1.2 1.2]) subplot(2,1,2); T = 0.1;f = 13; n = (0:T:1)'; xs = cos(2*pi*f*n); t = linspace(-0.5,1.5,500)'; ya = sinc((1/T)*t(:,ones(size(n))) - (1/T)*n(:,ones(size(t)))')*xs; plot(n,xs,'o',t,ya);grid; xlabel('Time, msec');ylabel('Amplitude'); title('Reconstructed continuous-time signal y_{a}(t)'); axis([0 1 -1.2 1.2]);
Power Law transformation A power law is a special kind of mathematical relationship between two quantities. When the frequency of an event varies as a power of some attribute of that event (e.g. its size), the frequency is said to follow a power law. For instance, the number of cities having a certain population size is found to vary as a power of the size of the population, and hence follows a power law. There is evidence that the distributions of a wide variety of physical, biological, and man-made phenomena follow a power law, including the sizes of earthquakes craters on the moon and of solar flares, the foraging pattern of various species the sizes of activity patterns of neuronal population the frequencies of words in most languages, frequencies of family names, the sizes of power outages and wars, and many other quantities.
The power law transformation’s equation is;
Practice Example
clc; clear all; I = imread('4107.ppm'); R =I(:,:,1); % red plane G =I(:,:,2); % green plane B =I(:,:,3); % blue plane %power law transformation on red plane image_double=im2double(R); [r c]=size(image_double); for i=1:r for j=1:c R1(i,j)=(1*power(image_double(i,j),0.4)); end end %power law transformation on green plane image_double=im2double(G); [r c]=size(image_double); for i=1:r for j=1:c G1(i,j)=(1*power(image_double(i,j),0.4)); end end %power law transformation on blue plane image_double=im2double(B); [r c]=size(image_double); for i=1:r for j=1:c B1(i,j)=(1*power(image_double(i,j),0.4)); end end % catenating three planes after transformation rgbb=cat(3,R1,G1,B1); % displaying images figure imshow(I) title('Given image') figure imshow(rgbb) title('Image after transformation') figure subplot(2,3,1) imshow(R) %plot title('RED intensity image') subplot(2,3,2) imshow(G); %plot title('GREEN intensity image') subplot(2,3,3) imshow(B); %plot title('BLUE intensity image') subplot(2,3,4) imhist(R) %plot title('RED intensity histogram') subplot(2,3,5) imhist(G); %plot
red intensity
Green intensity
Blue intensity
red intensity histogram Green intensity
