burg.m

burg 的matlab程序

clear

% N = 512;
% IP = 12;
% 
% t = 0:(N-1);
% fan1 = 0; fan2 = 0;
% fs = 100; %%%
% f1 = 0.2*fs;
% f2 = 0.25*fs;
% 
% x = sin(2*pi*f1/fs*t + fan1) + sin(2*pi*f2/fs*t + fan2);
% y = awgn(x,6);
% %  figure,plot(x); 
% %  figure,plot(y);
% x = y;


 x = [1,2,3,4,5];   %%%%% 做题
 N = length(x);
 IP =  2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% burg递推法

ef(0+1,:) = x;
eb(0+1,:) = x;
rou(0+1) = sum(x.^2)/N;  %% p0

p = 0;     %%%%%% 书上是从p=1开始

while p < IP   
    
    p = p + 1;
    
    use1 = ef(p,(p+1):N);
    use2 = eb(p,p:(N-1));
    k(p) = -2*sum(use1.*use2)/sum(use1.^2 + use2.^2);
    
    a(p,p) = k(p);
    rou(p+1) = (1 - k(p)^2)*rou(p-1+1);
    
%     rou2(p) = sum(use1.^2 + use2.^2)/(2*(N-p));   %%%%% 做题
%     rou3(p) = sum(use1.^2 + use2.^2)/2;
    
    for i = 1 : p-1
        a(p,i) = a(p-1,i) + k(p)*a(p-1,p-i);        
    end
    
%     ef(p+1,p+1:N) = ef(p-1+1,p+1:N) + k(p)*eb(p-1+1,p:N-1);     %% p+1 : N
%     eb(p+1,p+1:N) = eb(p-1+1,p:N-1) + k(p)*ef(p-1+1,p+1:N);     %% p+1 : N
    
    ef(p+1,p+2:N) = ef(p-1+1,p+2:N) + k(p)*eb(p-1+1,p+1:N-1);     %% p+2 : N
    eb(p+1,p+1:N-1) = eb(p-1+1,p:N-2) + k(p)*ef(p-1+1,p+1:N-1);   %% p+1 : N-1
    
end

% p = p + 1;      %%%%% 做题
% use1 = ef(p,(p+1):N);
% use2 = eb(p,p:(N-1));
% rou2(p) = sum(use1.^2 + use2.^2)/(2*(N-p));
% rou3(p) = sum(use1.^2 + use2.^2)/2;
 
a    %% 显示系数
[ef,eb]   %% 显示前向、后向预测误差(中间结果)
test = [rou;rou2;(rou - rou2)*100./rou;rou3]   %% 对比中间结果 %%%%% 做题


% test(:,3:6)*100
% a = [0.3678, 0.4612];
% b = [0.1439, 0.2365];
% 
% sum(a.^2 + b.^2)/2


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% 画频谱图
% p = IP;
% 
% M = 512;  %% 计算M点频谱
% t = 0:(M-1);
% 
% w = 2*pi*t/M;
% 
% for n = 1:M
%     temp = 0;
%     for i = 1:p
%         temp = temp + a(p,i)*exp(-w(n)*i*j);
%     end
%     yy(n) = rou(p+1)/abs(1+temp)^2;
% end
% 
% figure,plot([0:M-1]*2*pi/M - pi,fftshift(10*log10(yy)));
% title('burg法(dB)');
% 
% 
% % figure,plot([0:M-1]*20/M - 10,fftshift(10*log10(yy)));
% % grid on;
% % title('burg法(dB)');