democorr.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


%DEMO ON CORRELATION
clc,subplot,hold off,echo on
%                       CONCEPTS IN CORRELATION

%Correlation is a measure of similarity between two signals x(t) and y(t).

%Correlation is similar to convolution except that one signal is folded.

%The correlation of x(t) and y(t) is just the convolution of x(t) and y(-t)


%AUTO-CORRELATION
%The correlation of x(t) with itself is called autocorrelation Rxx

%The autocorrelation is an EVEN function with a maximum at the origin

%Reason: Autocorrelation of x(t) and x(t) is convolution of x(t) and x(-t).

%The maximum overlap (and maximum convolution) always occurs at t=0.

pause, % strike a key to continue

%EXAMPLE: Consider the sequence {x[n]}={0 1 2 3...8} described in MATLAB by
x=0:8

%To find the correlation, we use the routione corrnum. Here is what it does
help corrnum
%If we invoke this routine and plot rxx, here is what we get
Nx=length(x)-1;n=-Nx:Nx;rxx=corrnum(x,x)

pause     %Strike a key for the plot
dtplot(n,rxx,'o'),title('rxx'),pause    %This routine plots discrete signals

%Note the even symmetry of rxx and the maximum at the origin

%Meaning: Perfect correlation when the two signals are aligned
%	  Decreasing correlation on either side

pause     %Strike a key to continue
clc


%CROSS-CORRELATION
%The correlation of x(t) with y(t) is called crosscorrelation.

%Rxy and Ryx are not the same but are related.

%To see the relation, consider the signal
y=0:3
%Let us generate rxy and ryx to see what we get
rxy=corrnum(x,y);
ryx=corrnum(y,x);
[rxy;ryx]
%Do you see the difference? The two are actually folded versions of each other.
%This means that rxy and ryx start at different instants.

%If both sequences start at n=0, then
%rxy is the convolution x[n]*y[-n] and starts at n=-3
%ryx is the convolution y[n]*x[-n] and starts at n=-5

%Here is how we plot the two correlation sequences
Nx=length(x)-1;Ny=length(y)-1;nxy=-Nx:Ny;nyx=-Ny:Nx;mx=max([Nx Ny]);
my=1.2*max([rxy,ryx]);ax=[-mx mx 0 my];

pause     %Strike a key for the plot
subplot(211),axis(ax),dtplot(nxy,rxy,'o'),title('rxy'),...
subplot(212),axis(ax),dtplot(nyx,ryx,'o'),title('ryx'),pause

%Thus, we see that rxy[n]=ryx[-n]. Also, both have a maximum at n=0.

subplot,pause  % strike a key to continue
clc

%ANOTHER EXAMPLE OF CROSS CORRELATION
%If a signal x is correlated with its shifted version what do we get?
%Consider the signal x[n] and its shifted version x[n-4]
d=5;x=0:10,y=[zeros(1,d) x]  %y[n]=x[n-5]
%Here is what we get for the correlations
rxy=corrnum(x,y);ryx=corrnum(y,x);

Nx=length(x)-1;Ny=length(y)-1;nxy=-Nx-d:Ny-d;nyx=-Ny+d:Nx+d;
mx=max([Nx+d Ny+d]);my=1.2*max([rxy,ryx]);ax=[-mx mx 0 my];

pause     %Strike a key for the plot
subplot(211),axis(ax),dtplot(nxy,rxy,'o'),title('rxy'),...
subplot(212),axis(ax),dtplot(nyx,ryx,'o'),title('ryx'),pause
subplot,hold off

%NOTE: that rxy or ryx are just shifted versions of the autocorrelation rxx

%The shift (delay) is indicated by the location of the peak in rxy!

pause  % strike a key to continue


%APPLICATION: TARGET RANGING

%The above result is the basis for target ranging by radar (or ultrasound)

%We send a signal x which is reflected by the target (an airplane say)

%Ideally, the received signal is just a delayed version x(t-t0)

%The delay t0= 2Rc where R=target distance and c =velocity (of light)

%We crosscorrelate x(t) and x(t-t0) and locate the peak.

%This location gives t0. From t0, we estimate the target range R
 
pause   %strike a key to continue
clc
%OTHER APPLICATIONS OF CORRELATION

%Correlation is an effective method of finding signals buried in noise
%Noise is essentially uncorrelated. This means that if we correlate a noisy
%signal with itself, the correlation will be due only to any signal.
%This will exhibit a maximum at n=0 with decay on either side.

%The spectrum of the correlated noisy signal is much more useful than
%the spectrum of the noisy signal itself.

%Here are some examples
%Correlation of a noisy sine wave
%We generate a 20Hz sinusoid x for 4s, add noise n1 and nn2. The plots show:
%1. portions of the actual and noisy signals
%2. The various correlations
%3. spectra of the various correlations and the spectrum of the noisy signal 

%NOTE: Strike a key between plots. This may take a while
echo off
clear
ts=0.01;T=2;t=0:ts:T;N=length(t)-1;x=sin(2*20*pi*t);[zm,zn]=size(t);
if exist('version')==5,z1=eval('2*randn(zm,zn)');z2=eval('2*randn(zm,zn)');
else,uu='normal';eval('rand(uu)');z1=2*rand(zm,zn);z2=2*rand(zm,zn);end
x1=x+z1;x2=x+z2;xs=[x;x1;x2];x0=[z1;z2];
cn=corrnum(z1,z2,ts);cc=corrnum(x1,x2,ts);ca=corrnum(x1,x1,ts);
n=1:40;t=t-T/2;t1=-T:ts:T;f=(0:2*N)/(2*T);f1=(0:N)/T;
nspec=abs(fft(cn))/N;ccspec=abs(fft(cc))/N;
acspec=abs(fft(ca))/N;sfft=2*abs(fft(x1))/N;

subplot(221),plot(t(n),x(n)),title('noisefree signal x'),...
subplot(223),plot(t(n),z1(n)),title('noise signal n1'),...
subplot(222),plot(t(n),x1(n)),title('signal xn1 = x + n1'),...
subplot(224),plot(t(n),x2(n)),title('signal xn2 = x + n2'),pause

vx=matverch;
if vx < 4, eval('clg');else,eval('clf');end
subplot(221),plot(t1,cn,'-'),title('Correlation of n1 & n2'),...
subplot(222),plot(t1,ca,'-'),title('Autocorrelation of xn1 & xn1'),...
subplot(224),plot(t1,cc,'-'),title('Crosscorrelation of xn1 & xn2'),pause
subplot

subplot(221),plot(f(1:N+1),nspec(1:N+1)),...
title('Spectrum of noise correlation'),subplot(224),...
plot(f(1:N+1),ccspec(1:N+1)),title('Spectrum of crosscorrelation'),...
subplot(222),plot(f(1:N+1),acspec(1:N+1)),...
title('Spectrum of autocorrelation'),subplot(223),...
plot(f1(1:N/2+1),sfft(1:N/2+1)),title('FFT Spectrum of signal'),pause
subplot
echo on

%Note that the crosscorrelation and autocorrelation are a much better
%indicator of the actual signal than is the spectrum of the signal itself!

pause   % strike a key to continue
clc

%ANOTHER EXAMPLE
%Sinusoids buried in noise.

%Once again the crosscorrelation and autocorrelation are a much better
%indicator of the actual signal than is the spectrum of the signal itself!

%This is very important especially in situations where the signal strength 
%is not high or where the Signal to Noise Ratio (SNR) is very low.

%Strike a key between plots. This may take a while!!

echo off
clear
ts=.002;vr=version;if vr(1)=='S',ts=0.004;end
T=1;t=0:ts:T;N=length(t)-1;[zm,zn]=size(t);
x=2+2.5*sin(2*pi*90*t)-1.4*cos(2*pi*120*t)-4.5*sin(2*pi*220*t);
if exist('version')==5,z1=eval('4*randn(zm,zn)');
z2=eval('4*randn(zm,zn)');else,uu='normal';ru='uniform';
eval('rand(uu)');z1=4*rand(zm,zn);z2=4*rand(zm,zn);eval('rand(ru)');end
x1=x+z1;x2=x+z2;xs=[x;x1;x2];x0=[z1;z2];
cnn=corrnum(z1,z2,ts);css=corrnum(x1,x2,ts);
csn=corrnum(x1,z2,ts);ca1=corrnum(x1,x1,ts);
t1=-T:ts:T;scnn=abs(fft(cnn))/N;scss=abs(fft(css))/N;
scsn=abs(fft(csn))/N;sca1=abs(fft(ca1))/N;
f=(0:2*N)/(2*T);sfft=2*abs(fft(x1))/N;f1=(0:N)/T;
subplot(211),plot(t1,cnn,'-'),title('Correlation of two noise signals'),...
subplot(212),plot(t1,css,'-'),...
title('Crosscorrelation of two noisy signals'),pause

vx=matverch;
if vx < 4, eval('clg');else,eval('clf');end

subplot(221),plot(f(1:N+1),scnn(1:N+1)),title('Spectrum of noise n1*n2'),...
subplot(222),plot(f(1:N+1),sca1(1:N+1)),title('Spectrum of ns1*ns1'),...
subplot(224),plot(f(1:N+1),scss(1:N+1)),title('Spectrum of ns1*ns2'),...
subplot(223),plot(f1(1:N/2+1),sfft(1:N/2+1)),...
title('Spectrum of noisy signal'),pause
subplot
echo on

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