fig9_20.m

数字信号处理Matlab演示文件,其中各个文件加放置了不同的matlab子文件

%	Figure 9.20: Comparison of spectral estimates
%	Requires M-file arspe.m
%
%	BM
%	Jun 1998	


clear
clf
colordef(1,'black')

M = 64;				%no. of data points
num = 20;			%no, of realisations

%	sinewaves in noise
f = [ 8.5 10.5 12.5 ]/M;	%frequencies
a0 = [ -10 0 -10 ];		%relative power (dB)
phi = rand(size(f))*2*pi;	%random phase
npwr = -40;			%noise power in dB
t = 0:M-1;
a = sqrt(2)* 10 .^ (a0/20);
npwr = 10 ^ (npwr/20);
x = a*sin(2*pi*f'*t + diag(phi)*ones(size(f'*t)) )+npwr*randn(size(t));


%	set up data plot
subplot(211);
lo = line('Xdata',t,'Ydata',x,'Linestyle','o','Color','y','Erasemode','Xor');
ll = zeros(1,M);
for ii = 1:M
	ll(ii) = line('Xdata',[ii-1 ii-1],'Ydata',[0 x(ii)],'Linestyle','-','Color','y','Erasemode','Xor'); 
end
lh = line('Xdata',[0 M],'Ydata',[ 0 0 ],'Linestyle','-','Color','w','Erasemode','Xor');
xlabel('\itn'),ylabel('\itx(n)'),title('signal')
axis( [ 0 M -4 4 ]);

%	set up PSD plots
subplot(212)

%	DFT
N = 512;
H = hamming(M);
H = H/sqrt(H'*H);
xw = H'.* x;
Nfft = N;
X = fft(xw,Nfft);
ww = 0:Nfft/2;
ww = ww/Nfft;
psd1 = 10*log10(X.*conj(X));
ldft = line('Xdata',ww,'Ydata',psd1(1:Nfft/2+1),'Linestyle',':','Color','y','Erasemode','Xor');

%	AR covariance form
P = 8;
[ar,mmse,psd2] = arspe( x, P, Nfft );
lar = line('Xdata',ww,'Ydata',psd2(1:Nfft/2+1),'Linestyle','-','Color','m','Erasemode','Xor');
xlabel('normalised frequency'),ylabel('PSD (dB)')
axis( [ 0 0.5 -60 30 ])

%	position of sinewaves
lf1 = line( 'Xdata',[ f(1) f(1) ],'Ydata', [ -60 a0(1) ],'Linestyle','-','Color','w','Erasemode','Xor' );
lf2 = line( 'Xdata',[ f(2) f(2) ],'Ydata',[ -60 a0(2) ],'Linestyle','-','Color','w','Erasemode','Xor' );
lf3 = line( 'Xdata',[ f(3) f(3) ],'Ydata', [ -60 a0(3) ],'Linestyle','-','Color','w','Erasemode','Xor' );

fprintf(1,'Figure 9.20: press return to start\n')
pause


%	num realisation of signal
for ii=1:num

%	timer for graphics
t0 = clock;
te = 0;

%	refresh data
	phi = rand(size(f))*2*pi;	
	x = a*sin(2*pi*f'*t + diag(phi)*ones(size(f'*t)) )+npwr*randn(size(t));
	set(lo,'Xdata',t,'Ydata',x);
	for jj=1:M
		set(ll(jj),'Xdata',[jj-1 jj-1],'Ydata',[0 x(jj)]);
	end
	
%	update DFT
	X = fft((H' .* x),Nfft);
	psd1 = 10*log10(X.*conj(X));
	set(ldft,'Xdata',ww,'Ydata',psd1(1:Nfft/2+1));

%	update AR spectral estimate
	[ar,mmse,psd2] = arspe( x, P, Nfft );
	set(lar,'Xdata',ww,'Ydata',psd2(1:Nfft/2+1));

%	constant frame rate
while te < 2.0
	te = (clock - t0) * [0 0 86400 3600 60 1]';
end

end