hpb1.m

谱估计及阵列信号处理算法仿真库

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%   (a)    %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Generate the autocorrelation function
M = 10;
r = zeros(M+1,1);
for k=1:M+1,
	r(k) = (M-(k-1))/M;
end

rr = [r(M+1:-1:1); r(2:M+1)];
figure(1);
plot(rr)

%%%%%%%%%%%%%%%%
%Compute the PSD
L = 256;
x = [r(1:M+1); zeros(L-(2*M+1),1); r(M+1:-1:2)];
size(x)

xf = fft(x);
xf_real = real(xf);
xf_imag = imag(xf);

figure(2); plot(xf_real);
figure(3); plot(xf_imag);
xf = xf_real;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The effect of truncation of imaginary part
x1 = [r(M+1:-1:2); r(1:M+1); zeros(L-(2*M+1),1)];
size(x1)

xf1 = fft(x1);
xf_real1 = real(xf1);
xf_imag1 = imag(xf1);

figure(4); plot(xf_real1);
figure(5); plot(xf_imag1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%   (b)    %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Analytic results for the PSD
yf = zeros(L,1);
yf(1) = M;
for f=2:L,
	omega = (f-1)*2*pi/L;
	yf(f) = 1/M*(sin(omega*M/2)/sin(omega/2))^2;
end
y = real(ifft(yf));
figure(6)
plot(1:L,y,'b:',1:L,x,'r');
plot(y)
figure(7)
plot(abs(x-y))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%   (c)    %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%
%Yule-Walker equations
a1 = -0.9;
rr = [1 a1; a1 1]\[1; 0]
r = zeros(M+1,1);
r(1:2) = rr(1:2);
for k=3:M+1,
	r(k) = (-a1)^(k-2)*r(2);
end

rr = [r(M+1:-1:1); r(2:M+1)];
figure;
plot(rr)

%%%%%%%%%%%%%%%%%%%
%Residue method
r1 = zeros(M+1,1);
for k=1:M+1,
	r1(k) = (-a1)^(k-1)/(1-a1^2);
end
rr1 = [r1(M+1:-1:1); r1(2:M+1)];

figure;
plot(abs(rr-rr1))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The analytic values for the autocorrelation function
M = 10;
%M = 127;

r = zeros(L/2+1,1);
for k=1:L/2+1,
	r(k) = (-a1)^(k-1)/(1-a1^2);
end
% clear x
% x = [r(L/2+1:-1:1); r(2:L/2)];

x = [r(1:M+1); zeros(L-(2*M+1),1); r(M+1:-1:2)];
size(x)

xf = real(fft(x));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Analytic values for the PSD
L=256;
yf = zeros(L,1);
for f=1:L,
	omega = (f-1)*2*pi/L;
	yf(f) = 1/(1+a1^2+2*a1*cos(omega));
end

%Compare computed and analytic PSD
figure; plot(1:L,xf,'b',1:L,yf,'r')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Compute the autocorrelation function by usig ifft
y = real(ifft(yf));

figure; plot(1:L,x,'b',1:L,y,'r')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%   (d)    %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure; 
xf = fftshift(real(fft(x)));
plot(xf)

omega0 = 2*pi/L;
figure; plot(-pi:omega0:pi-omega0,xf)