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


vx=matverch;

clc,subplot,hold off,echo on
%                     SAMPLING  ALIASING  AND  LEAKAGE

%SAMPLING:
%To sample a continuous signal x(t) without loss of information, fs > 2*f0.
%fs is the sampling frequency
%f0 is the highest frequency in x(t). 

%The frequency 2*f0 is called the Nyquist rate or Nyquist frequency
%The sampling interval Ts=1/2f0 is called the Nyquist interval.
%To sample a pure sinusoid, we must acquire MORE THAN 2 samples/period.

%ALIASING:
%Undersampling results in the ALIASING of high frequencies to lower ones

%LEAKAGE
%If we sample periodic signals for non-integer periods, the spectral
%frequencies leak to neighbouring values. This is leakage.
%Leakage occurs even if we satisfy the Nyquist sampling rate.

pause   %strike a key to continue
clc
%PART 1: NO ALIASING

%Let  x(t) = 1.3*sin(200*pi*t)

%It is a pure sinusoid at 100Hz with a peak value of 1.3.
%Its twosided spectrum has magnitudes of 0.65 at +100Hz and -100Hz

%To avoid aliasing,we must take at least 2 samples per time period.
%Choose a sampling rate of 8 samples per 3 periods.

%Thus ts=3T/8=0.375T   and   fs = 266.66Hz


%Compute the spectrum of the sampled signal
f0=100;T=1/f0;		%Frequency and period		
ts=0.375*T;        	%Sampling interval [sampling rate =1/ts = 266.66hz]
n=32; 			%Choose 32 samples (i.e. sample x(t) 12 periods)
t=0:ts:(n-1)*ts;	%Time array 
df= 1/(ts*n);		%The frequency spacing be df= 1/(n*ts)
sf=[-16:15]*df; 	%Frequency array corresponding to ts
x=1.3*sin(2*f0*pi*t); 	%The sampled function
X=abs(fft(x))/32;  	%And its FFT magnitude
X=fftshift(X);   	%Re-order to give two sided spectrum
t1=0:0.1*ts:(n-1)*ts;x1=1.3*sin(2*f0*pi*t1);

pause  %strike a key to see the plots

subplot(211),dtplot(t,x,'o'),hold on,plot(t1,x1,'--'),grid,...
title('sampled time function'),hold off,subplot(212),dtplot(sf,X,'o'),...
title(['f0 = ',num2str(f0),'  sf = ',num2str(1/ts)]),pause

pause  %Strike a key to continue
clc
%PART 2: EFFECTS OF ALIASING

%x(t) = 1.3*sin(200*pi*t)    f0=100Hz

%To see aliasing, choose a rate LESS than 2 samples per period
%Take 4 samples per 3 periods. 

%Then ts=0.75T and  fs=133.33Hz

%Since f0=100Hz, the aliased component should appear at 133.33-100 = 33.33Hz

%Compute the spectrum of the sampled signal
ts=0.75*T;         	%Sampling rate = 1/ts = 133.33hz
t=0:ts:(n-1)*ts;      	%Time array to generate 32 samples
df=1/(ts*n);sf=[-16:15]*df; %Frequency array corresponding to ts
x=1.3*sin(2*f0*pi*t);X=fft(x)/32;  %Sampled function & its FFT
X=abs(fftshift(X));%Re-order magnitude spectrum
t1=0:0.1*ts:(n-1)*ts;x1=1.3*sin(2*f0*pi*t1);

pause   %strike a key to see the plots

if vx < 4, eval('clg');else,eval('clf');end
subplot(211),dtplot(t,x,'o'),hold on,plot(t1,x1,'--'),grid,...
title('sampled time function'),hold off,subplot(212),dtplot(sf,X,'o'),...
title(['f0 = ',num2str(f0),'  sf = ',num2str(1/ts)]),pause

pause    %strike a key to continue
clc

%PART 3: MORE ON ALIASING

%TWO SINUSOIDS at 100hz and 50hz 

%x(t) = 1.3*sin(200*pi*t) + 2.1*sin(100*pi*t)

%We choose fs = 133.33Hz 

%Only the 100hz component should be aliased (to 33.33hz)
%The 50Hz sinusoid (unaliased) should show up at 50Hz


x=1.3*sin(200*pi*t)+2.1*sin(100*pi*t);X=fft(x)/32;X=abs(fftshift(X));
x1=1.3*sin(200*pi*t1)+2.1*sin(100*pi*t1);

pause     %strike a key to see the plots

if vx < 4, eval('clg');else,eval('clf');end
subplot(211),dtplot(t,x,'o'),hold on,plot(t1,x1,'--'),grid,...
title('sampled time function'),hold off,subplot(212),dtplot(sf,X,'o'),...
title(['f1 = 50   f2 = 100   sf = ',num2str(1/ts)]),pause

pause    %strike a key to continue
clc

%TWO SINUSOIDS at 150hz and 225hz    fs=133.33Hz

%x(t) = 1.3*sin(300*pi*t) + 2.1*sin(450*pi*t);

%Both components should show aliasing.

%The 150Hz component will alias to (|150-133.33|)=16.66hz

%The 225Hz component will alias to |225-133.33K|=|225-266.67|=41.67Hz


x=1.3*sin(3*f0*pi*t)+2.1*sin(4.5*f0*pi*t);X=fft(x)/32;X=abs(fftshift(X));
x1=1.3*sin(3*f0*pi*t1)+2.1*sin(4.5*f0*pi*t1);

pause   %strike a key to see the plots

if vx < 4, eval('clg');else,eval('clf');end
subplot(211),dtplot(t,x,'o'),hold on,plot(t1,x1,'--'),grid,...
title('sampled time function'),hold off,subplot(212),dtplot(sf,X,'o'),...
title(['f1 = 150   f2 = 225   sf = ',num2str(1/ts)]),pause

pause    %strike a key to continue

clc

%THREE SINUSOIDS at 50hz and 100hz & 350hz    fs=133.33Hz

%x(t) = sin(100*pi*t) + 2sin(200*pi*t) - 3sin(700*pi*t)

%The 50Hz component will show no aliasing
%The 100Hz component will alias to 33.33Hz
%The 350hz component will alias to 50hz & ADD to the 50hz component!
%THE MAGNITUDE OF THE 50Hz COMPONENT SHOULD THUS EQUAL 2


x=sin(100*pi*t)+2*sin(200*pi*t)-3*sin(700*pi*t);X=fft(x)/32;X=abs(fftshift(X));
x1=sin(100*pi*t1)+2*sin(200*pi*t1)-3*sin(700*pi*t1);

pause    %strike a key to see the plots

if vx < 4, eval('clg');else,eval('clf');end
subplot(211),dtplot(t,x,'o'),hold on,plot(t1,x1,'--'),grid,...
title('sampled time function'),hold off,subplot(212),dtplot(sf,X,'o'),...
title(['f1=50  f2=100  f3=350  sf = ',num2str(1/ts)]),pause

pause   %strike a key to continue

clc

%PART 5: ALIASED SPECTRUM OF COMBINATIONS

%Consider a combination of cosines at multiples of 10Hz with

%f = 10 40 200 220 260 280 300 340 360 380 Hz  Magnitude = 1, phase = 0
%The period T=0.1s     Choose fs=160Hz (16 samples/period)

%The actual and aliased frequencies are
%10 40 200 220 260 280 300 340 360 380 Hz, Original fequencies 
%10 40  40  60  60  40  20  20  40  60 Hz, Aliased  frequencies

%The components at 20Hz, 40Hz  and 60Hz add up to give the sampled signal
%xs(t)= cos(2pi*10t) + 2cos(2pi*20t) + 4cos(2pi*40t) + 3cos(2pi*60t)

%Its spectrum has magnitudes :  0.5  1   2  and 1.5
%at the frequencies          :  10   20  40 and  60 hz
                
%We plot the spectra of the actual and sampled functions

pause    %strike a key to see the plots

k=[10 40 200 220 260 280 300 340 360 380]; %frequencies of components
t=0.1*(0:127)/128; %x=sum(cos(k'*2*pi*t));
x=0*t;for j=1:length(k);x=x+cos(k(j)*2*pi*t);end
ts=1/160;n=16;t1=0:ts:(n-1)*ts; %x1=sum(cos(k'*2*pi*t1));
x1=0*t1;for j=1:length(k);x1=x1+cos(k(j)*2*pi*t1);end 
df= 1/(ts*n);sf=[-8:7]*df;sf1=(-64:63)*10;
y=abs(fft(x))/128;y1=abs(fft(x1))/n;s=[-400 400 0 2];hold off
y=fftshift(y);y1=fftshift(y1);A='Actual Function';S='Sampled Function';
if vx < 4, eval('clg');else,eval('clf');end

subplot(221),plot(t,x),grid,title(A),subplot(222),plot(t,x),...
hold on,dtplot(t1,x1,'o'),grid,title(S),hold off,subplot(223),...
axis(s),dtplot(sf1,y),axis(s);title(['spectrum: ' A]),subplot(224),axis(s),...
dtplot(sf,y1),axis(s),title(['spectrum: ' S]),pause


pause    %strike a key to continue

clc

%PART 6

%LEAKAGE EFFECTS

% Let x(t)=sin(5*pi*t)

%f0=2.5Hz,  Period T=0.4s

%We sample for three signal durations: T1=T   T2=4.5T   T3=21.5T

%We should see no leakage for T1 (one full period)

%The leakage effect for T3 (longer duration) should be less than that for T2

%Choose fs=10Hz

pause    %Strike a key for plots

x='sin(5*pi*t)';
td1 = 0.4000;
td2 = 1.8000;
td3 = 8.6000;
t=0:.1:td1-.1;x1=eval(x);
t=0:.1:td2-.1;x2=eval(x);
t=0:.1:td3-.1;x3=eval(x);
y1=abs(fftshift(fft(x1)))/4;y2=abs(fftshift(fft(x2)))/18;
y3=abs(fftshift(fft(x3)))/86;
f1=(-2:1)/td1;f2=(-9:8)/td2;f3=(-43:42)/td3;
A='o';B='--';D='.';E='-';G='Duration = ';H='o';
subplot,hold off,subplot(221),s=[-5 5 0 .6];
dtplot(f1,y1,H),axis(s);title([G 'T']),subplot(222),axis(s),dtplot(f2,y2,A),...
axis(s);title([G '4.5T']),subplot(223),axis(s),dtplot(f3,y3,D),axis(s);...
title([G '21.5T']),subplot(224),axis(s),dtplot(f1,y1),axis(s);...
hold on,dtplot(f2,y2,A),plot(f2,y2,'-'),dtplot(f3,y3,D),plot(f3,y3,'-'),...
title('Overplot of Results'),pause

%We see no leakage for T1 (one full period)
%The leakage effect for T3 (longer duration) is less than that for T2

pause  %Strike a key to continue

hold off,subplot,clc

%			ALIASING & LEAKAGE: A MOVIE

%We show how the aliasing component varies as the frequency of a
%sinusoid is changed from 225hz to 450hz keeping fs=120Hz
%Since we have not sampled for integer periods, Leakage is also present

%Note how -0.5fs & 0.5fs (-60Hz & 60Hz) act as a wall
%for the highest frequency we can observe in the spectrum

pause %Strike any key to see the movie
echo off
clear
vx=matverch;if vx>=4,ml=4;er='erasemode';bk='back';xd='xdata';yd='ydata';
else,ml=3;end
subplot,hold off
if vx < 4, eval('clg');else,eval('clf');end

s=[-60 60 0 1];
if ml==4,plot([0 0.5],[1 1],'--'),axis(s);hold on,eval('grid on'),else
axis(s);plot([-60],[1],'.i',60,1,'.i'),grid,hold on,end  %%%NEW

X1=0;sf1=0;if ml==4,plt=eval('plot(sf1,X1,er,bk)');end
f0=100;ts=1/120;nn=64;t=0:ts:(nn-1)*ts;df=1/(ts*nn);sf=[-nn/2:nn/2-1]*df;
if ml==4,dk=0.01;else,dk=0.05;end
for k=4.5:dk:9,
x=1.3*sin(k*f0*pi*t);X=fft(x)/nn;X=abs(fftshift(X));
if ml==4,eval('set(plt,xd,sf,yd,X);drawnow'),else
plot(sf1,X1,'i'),plot(sf,X),grid,sf1=sf;X1=X;end
end,if ml==3,pause,end,hold off
echo on

pause  %STRIKE A KEY TO END THIS DEMO
echo off
clear,clc