基于MATLAB的谱估计和自适应程序 求助
 

以下为谱估计和自适应程序，哪位大神能帮我注解下程序的语句和告诉我流程图是什么？
clear;

clc;

%Burg Data leng N=64

clear;

clc;

N=1024;

n=[0:N-1];

randn('state',0);

wn=randn(1,N);

xn=sqrt(20)*sin(2*pi*0.1*n)+sqrt(20)*sin(2*pi*0.4*n)+wn;

ef(1,:)=xn;

eb(1,:)=xn;

tempsum=0;

for i=1:1:N

tempsum=tempsum+xn(i).*xn(i);

end

thegma(1)=tempsum/N;



%p=1

sum1=0;

sum2=0;

for n=2:N,

sum1=sum1+ef(1,n)*eb(1,n-1);

sum2=sum2+abs(ef(1,n))^2+abs(eb(1,n-1))^2;

end

a(1,1)=-2*sum1/sum2;

thegma(2)=(1-abs(a(1,1))^2)*thegma(1);

for n=3:N,

ef(2,n)=ef(1,n)+a(1,1)*eb(1,n-1);

eb(2,n)=eb(1,n-1)+a(1,1)*ef(1,n);

end



p=1;

IP=14;

while p<IP

p=p+1;

sum1=0;

sum2=0;

for i=p+1:N

sum1=sum1+ ef(p,i)*eb(p,i-1);

sum2=sum2+ef(p,i)*ef(p,i)+eb(p,i-1)*eb(p,i-1);

a(p,p)=-2*sum1/sum2;

end

thegma(p+1)=[1-(abs(a(p,p)))^2]*thegma(p);

for i=p+2:N

ef(p+1,i)=ef(p,i)+a(p,p)*eb(p,i-1);

eb(p+1,i)=eb(p,i-1)+a(p,p)*ef(p,i);

end

for k=1:1:p-1

a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);

end

end

b=[1,a(p,:),zeros(1,N-p-1)];

Pxx=abs(thegma(p))./((abs(fft(b))).^2);

Pxx=Pxx(1:N/2);

Pxx=10*log10(Pxx);

n=[0:N-1];

n=n(1:N/2);

% subplot(2,2,1);

plot(n/N,Pxx);

xlabel('Frequency');

ylabel('Power Spectrum (dB)');

title('Burg Data leng N=64');

grid on;




2、自适应的程序

clear all

close all

%channel system order

sysorder = 5 ;

% Number of system points

N=2000;

inp = randn(N,1);

n = randn(N,1);

[b,a] = butter(2,0.25);

Gz = tf(b,a,-1);

%This function is submitted to make inverse Z-transform (Matlab central file exchange)

%The first sysorder weight value

%h=ldiv(b,a,sysorder)';

% if you use ldiv this will give h :filter weights to be

h= [0.0976;

0.2873;

0.3360;

0.2210;

0.0964;];

y = lsim(Gz,inp);

%add some noise

n = n * std(y)/(10*std(n));

d = y + 3*n;

totallength=size(d,1);

%Take 60 points for training

N=60 ;

%begin of algorithm

w = zeros ( sysorder , 1 ) ;

for n = sysorder : N

u = inp(n:-1:n-sysorder+1) ;

y(n)= w' * u;

e(n) = d(n) - y(n) ;



if n < 20

mu=0.32;

else

mu=0.15;

end

w = w + mu * u * e(n) ;

end





%check of results

for n = N+1 : totallength

u = inp(n:-1:n-sysorder+1) ;

y(n) = w' * u ;

e(n) = d(n) - y(n) ;

end

hold on

plot(d)

plot(y,'r');

title('System output') ;

xlabel('Samples')

ylabel('True and estimated output')

figure

semilogy((abs(e))) ;

title('Error curve') ;

xlabel('Samples')

ylabel('Error value')

figure

plot(h, 'k+')

hold on

plot(w, 'r*')

legend('Actual weights','Estimated weights')

title('Comparison of the actual weights and the estimated weights') ;

axis([0 6 0.05 0.35])



% RLS 算法

randn('seed', 0) ;

rand('seed', 0) ;



NoOfData = 8000 ;

Order = 32 ;

Lambda = 0.98 ;

Delta = 0.001 ;

x = randn(NoOfData, 1) ;

h = rand(Order, 1) ;

d = filter(h, 1, x) ;



P = Delta * eye ( Order, Order ) ;

w = zeros ( Order, 1 ) ;



% RLS Adaptation



for n = Order : NoOfData ;



u = x(n:-1:n-Order+1) ;

pi_ = u' * P ;

k = Lambda + pi_ * u ;

K = pi_'/k;

e(n) = d(n) - w' * u ;

w = w + K * e(n) ;

PPrime = K * pi_ ;

P = ( P - PPrime ) / Lambda ;

w_err(n) = norm(h - w) ;



end ;



% Plot results



figure ;

plot(20*log10(abs(e))) ;

title('Learning Curve') ;

xlabel('Iteration Number') ;

ylabel('Output Estimation Error in dB') ;



figure ;

semilogy(w_err) ;

title('Weight Estimation Error') ;

xlabel('Iteration Number') ;

ylabel('Weight Error in dB') ;