democonv.m

很多matlab的源代码

% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


clc,subplot
vx=matverch;
% if vx < 4, eval('clg');else,eval('clf');end
echo on
%CONVOLUTION USING MATLAB
%Being a numerical program, MATLAB can only approximate the convolution of 
%continuous time signals. The smaller the interval at which the two signals 
%are described, the more accurate the results.

%EXAMPLE:
%Let us convolve a cosine pulse with a triangular pulse 
%Describe the signals x and y using the following MATLAB statements:
ts=0.01;     		%Choose an interval.
t1=-1:ts:1;           	%Choose a time axis for x.
x=3*cos(pi*t1/2);     	%Choose the signal x as a cosine pulse.
t2=0:ts:2;           	%Choose a time axis for y.
y=2*uramp(t2)-4*uramp(t2-1)+2*uramp(t2-2);  %signal y is a triangular pulse.
ax=[-1 3 0 5];

pause  %Strike a key to see the signals x and y
axis(ax);
plot(t1,x,t2,y),axis(ax),grid,title('The signals x and y'),pause
%if v==3,pause,end
hold on

pause  %strike a key to continue 

%The convolution z=x*y is obtained using the command:
z1=conv(x,y);         %Find the convolution z of x and y.

%MATLAB does not use the interval ts. Nor does it tell us the correct start. 
%The MATLAB result must therefore be modified as follows:
%1.   Include the interval. The true convolution z equals ts*z1:
%2.   Establish the correct range. The convolution extends from -1 to 3.
%To plot z, we need a time axis between -1 and 6 at intervals of  ts.
%These modifications are described by the statements:

z=ts*z1;                 %Modify the result by including the interval.
t3=-1:ts:3;              %Create a time axis for z.

pause %strike a key for the convolution plot
plot(t3,z),title('x and y and z=x*y'),pause
echo off,hold off
if vx>=4,xx='2*uramp(t)-4*uramp(t-1)+2*uramp(t-2)';txx=[0 2];
yy='3*cos(pi*t/2)';tyy=[-1 1];eval('convplot(xx,yy,txx,tyy)');end
echo on
pause  %strike a key to continue 

%SOME PROPERTIES OF CONVOLUTION USING MATLAB
%As a consistency check, we can verify the area property as follows:
ax=sum(x)*ts;            %Area under x.
ay=sum(y)*ts;            %Area under y.
az=sum(z)*ts;            %Area under z.
axy=ax*ay;               %Product of areas under x and y.
df=az-axy;               %Their difference.
p=[ax;ay;az;axy;df]      %Create a vector of the results and display.

pause  %strike a key to continue 

%To finding the maximum convolution and its time of occurrence, use:
zmax=max(z)              %Display the maximum values of z
k=find(z==zmax);         %Find index corresponding to zmax; Note the ==
tmax=-1+ts*(k-1);        %Display tmax. Use k-1 (In Matlab index starts at 1)
disp([zmax tmax])        %and add -1 unit since the conv starts at -1.

pause  %strike a key to continue 

%CONVOLUTION OF SIGNALS WITH DISCONTINUITIES
%For signals with jumps, what value do we choose at the discontinuities?
%The midpoint is a reasonable choice.
%This means using the trapezoidal rule for computing areas (and convolution)
%Consider the convolution of two pulses of unit height and width.
%Their convolution is a triangle with width 2 and height 1 with endpoints 0

%Here is what Matlab gives us if we use 1 for the endpoints of each pulse
ts=0.1;t1=0:ts:1;		%time interval and array
x=ones(1,11); 			%pulse sampled at 0.1 so 11 samples
y=x;				%same signal
z=ts*conv(x,y);t2=0:ts:2;	%Convolution
za=uramp(t2)-2*uramp(t2-1)+uramp(t2-2);	%Actual theoretical result

%We overplot z and za to see any differences

pause  %strike a key to see the plot
plot(t2,z,t2,za),title('Actual/numerical convolution of two pulses'),pause

pause  %strike a key to continue 


%To improve the results, modify x and y by choosing 0.5 for the endpoints.
%Here is one way to do this
x=(t1>0 & t1<1)+.5*(t1==0 | t1==1);y=x;

%Recompute the convolution
z2=ts*conv(x,y);

pause  %strike a key to see the plot
plot(t2,z2,t2,za,'--',t2,z,':'),...
title('Conv. of pulses with (a)endpoints=1 (b)endpoints=0.5 (c)Actual'),pause
pause  %strike a key to continue 


%CHOOSING A SAMPLING INTERVAL
%What sampling interval should we choose for "good"  results? 
%We could use the sampling theorem (choose ts=1/2Fh where Fh=highest freq) 
%Problem: We do not know Fh (it is actually infinite for timelimited signals!)

%A PRACTICAL APPROACH
%Another approach is to choose ts to minimize some error criterion.
%A useful criterion is based on minimizing the RELATIVE energy error.
%Here is an iterative approach to establish the sampling interval and 
%number of samples n based on two signals with equal durations

%1.   	Start with a small n. Find the energy E1 in the convolution.
%2.   	Increase n to n+1 and find the energy E2 in the new convolution.
%3.   	Find the relative error D=abs(E2-E1)/E1.
%4.	Compare with the required relative error, say D0.
%5.   	Repeat Steps 2 - 4 until D is less than D0.

%Do you see the advantages?
%First, we do not need the actual energy in the convolution!!
%The sampling interval ts cancels out in the computation of D!!
%Finally ts = total signal duration/n

%The file converr.m shows how you can code this in Matlab
%It is set up to plot the energy error (5%) vs the number of samples for the 
%convolution of two identical rectangular pulses (with endpoints=1)

pause  %Strike a key to see how many samples are required for an error of < 5%
converr

pause  %strike a key to continue 
clc
%DISCRETE CONVOLUTION
%Discrete convolution uses the same statements as described previously.
%We choose an integer index instead of a time index. 
%For sampled convolution, the result must be multiplied by the interval ts.
%Here is an example where the sampling interval is unity:

n1=0:10;                %Choose a time index for x.
x=2*(.8).^(n1);         %Choose the signal x.
n2=-4:4;                %Choose a time index for y.
y=2*sin(pi*n2/4);       %Choose the signal y.
z=conv(x,y);            %Find the convolution z of x and y.
n3=-4:14;               %Create a time index for z.

pause  %strike a key to see the plot
subplot(221),dtplot(n1,x,'o'),axesn,title('signal x'),...
subplot(222),dtplot(n2,y,'o'),axesn,title('signal y'),...
subplot(223),dtplot(n3,z,'o'),axesn,title('Convolution x*y'),pause,subplot

pause  %strike a key to continue 

%PERIODIC CONVOLUTION
%MATLAB provides no routines for computing periodic convolution.  We  supply 
%the routines convp for periodic convolution using 3 different approaches:

%1.   convp(x,y,dt) uses the wraparound method (Default)
%2.   convp(x,y,dt,'c') uses the circulant matrix.
%3.   convp(x,y,dt,'f') uses the FFT algorithm.

%The interval dt automatically scales  the  convolution  if  specified.  For 
%plotting results, the convolution range must be  established.  The  results 
%from each are identical to within roundoff error.
%We provide an example using the wraparound routine assuming both  sequences 
%to start at 0 with an interval of 0.5 between samples:

x=[2 1 3 6 2];           %Choose the signal x with 5 samples.
y=[3 2 1 2 3];           %Choose the signal y with 5 samples.
z=convp(x,y,0.5);       %Find the periodic convolution z of x and y.
n=0:.5:2;                %Time index for plotting z.

pause, %strike a key for plot
ax=[0 7 0 60];
 if vx < 4, eval('clg');else,eval('clf');end
subplot(221),dtplot(n,x,'o'),axesn,title('One period of signal x'),...
subplot(222),dtplot(n,y,'o'),axesn,title('One period of signal y'),...
subplot(212),dtplot(n,z,'o');axis(ax),axesn,...
title('Periodic conv of 1 period of x and y'),pause
%Plot the periodic convolution z.

%The other two approaches employ a similar syntax.

pause  %strike a key to continue 


%EFFECT OF DURATION ON PERIODIC CONVOLUTION
%The results of periodic convolution depend on how many  periods  we  choose 
%for the convolution. If we use three periods,  the  convolution  result  is 
%scaled threefold. Here is the previous example using three periods:

x3=[x x x];              %Choose three periods of signal x, 15 samples.
y3=[y y y];              %And three periods of signal y, 15 samples.
z3=convp(x3,y3,0.5);    %Find periodic convolution over three periods.
n1=0:.5:7;               %Time index for plotting z(0-2  2.5-4.5  5-7).

pause  %strike a key to see the plot
if vx < 4, eval('clg');else,eval('clf');end

subplot(211),dtplot(n1,z3,'o');axis(ax);axesn,...
title('Periodic convolution of 3 periods of x and y'),...
subplot(212),dtplot(n,z,'o');axis(ax),axesn,...
title('Periodic conv of 1 period of x and y'),pause
subplot

%Note that

%The first 5 elements in z3 are three times the convolution z 
%And  repeated three times.
%It is therefore common to normalize the convolution results by 
%the number of periods.

pause  %strike a key to continue 

%THE CENTRAL LIMIT THEOREM
%The repeated convolution of energy signals tends to a  Gaussian  shape.  We 
%supply the routine clt that demonstrate the central limit theorem. Type

%clt(ty,n);

%This gives the convolution of n signals of type ty
%ty='t' for tri, 'r' for rect or 'e' for exponentials.

%Each routine compares the convolution with the expected Gaussian result.
%For the repeated convolution of exponentials 20 times, here is the result: 
clt('e',20);pause

pause  %strike a key to continue 
clc

%CORRELATION
%The correlation Rxy of x and y is just the convolution of x and the  folded 
%version of y. To get Ryx, we perform the convolution of the folded  version 
%of x with the signal y. Remember that Rxy and Ryx are not the same!! 
%Here is how we could find the correlation of two signals:

x=[2 1 3 6];             %Choose the signal x with 4 samples.
y=[3 2 1 2 3];           %Choose the signal y with 5 samples.
xf=fliplr(x);            %Fold the signal x.
yf=fliplr(y);            %Fold the signal y.
rxy=conv(x,yf);          %Find the correlation Rxy of x and y.
ryx=conv(y,xf);          %Find the correlation Ryx of y and x.
n1=-3:4;                 %Time index for plotting Rxy.
n2=-4:3;                 %Time index for plotting Ryx.

pause  %strike a key to see the plot
subplot(221),dtplot(0:3,x,'o'),axesn,title('signal x'),...
subplot(222),dtplot(0:4,y,'o'),axesn,title('signal y'),...
subplot(223),dtplot(n1,rxy,'o'),axesn,title('correlation Rxy'),...
subplot(224),dtplot(n2,ryx,'o');axesn,title('correlation Ryx'),pause
subplot

%Note that Rxy and Ryx are folded versions of each other
[rxy;ryx]

%There is more about correlation in the m-file democorr.m

pause  %strike a key to END THIS DEMO 
subplot,hold off,echo off
clear