bt.m

there is a simple demo for non parameteric spectral estimation methods

function Px = bt(x)

global sigtype;
global vari;
global ar_coef;
global freq;
global amp;

HH = findobj(gcf,'Tag','bt_n');
p = str2num(get(HH,'String'));
if ~p
    error('Invalid N');
else
    N=p;
end;

HH = findobj(gcf,'Tag','bt_m');
p = str2num(get(HH,'String'));
if ~p
    error('Invalid M');
else
    M=p;
end;

HH = findobj(gcf,'Tag','bt_win');
winval = get(HH,'Value');

if(winval==2)
    win = 4;
elseif(winval==3)
    win = 3;
elseif(winval==4)
    win = 2;
elseif(winval==5)
    win = 5;
else
    win = 1;
end;

figure('Name','Blackman-Tukey Method','Menubar','none', ...
       'NumberTitle', 'off','Position',[212 200 450 500]);

TotN = 512;
j = sqrt(-1);
noise = sqrt(vari)*randn(50*N,1);
seta1 = (rand(1)-0.5)*2*pi;
seta2 = (rand(1)-0.5)*2*pi;

index_f = linspace(0,1,TotN);

Px_ave = 0;

subplot('Position',[0.13 0.581 0.77 0.344]);
for cnt=1:50; 
if(sigtype == 2)
    for nn = 1:N;
        sigdata(nn) = amp(1)*sin(nn*freq(1)*pi+seta1)+amp(2)*sin(nn*freq(2)*pi+seta2)+noise((cnt-1)*N+nn);
    end;
elseif(sigtype == 3)
    sigdata(1) = noise((cnt-1)*N+1);
    sigdata(2) = ar_coef(1)*sigdata(1)+noise((cnt-1)*N+2);
    sigdata(3) = ar_coef(1)*sigdata(2)+ar_coef(2)*sigdata(1)+noise((cnt-1)*N+3);
    for nn = 4:N;
        sigdata(nn) = ar_coef(1)*sigdata(nn-1)+ar_coef(2)*sigdata(nn-2)+ar_coef(3)*sigdata(nn-3)+noise((cnt-1)*N+nn);
    end;
else
    for nn = 1:N;
        sigdata(nn) = noise((cnt-1)*N+nn);
    end;
end;

newsig = sigdata - ones(1,N)*(sum(sigdata)/N);
fil_b = fliplr(newsig);
fil_out = (filter(fil_b,[1],newsig))/(N-1);
auto_sig = fliplr(fil_out);
R = auto_sig(2:M);
r  = [fliplr(R),auto_sig(1),R];
win_M  = 2*M-1;
w  = ones(win_M,1);
if (win == 2) w = hamming(win_M);
   elseif (win == 3) w = hanning(win_M);
   elseif (win == 4) w = bartlett(win_M);
   elseif (win == 5) w = blackman(win_M); 
   end;
r = r'.*w;
Px = abs(fft(r,1024));
Px(1)=Px(2);
Px_ave = Px_ave + Px/50.0;
Px_db = 10*log10(Px);
plot(index_f,Px_db(1:512));
hold on;
end;
hold off;

xlabel('Frequency (units of \pi)');
ylabel('Magnitude (dB)');
title('Overlay plot of 50 Blackman-Tukey estimates','fontsize',13);

Px_ave_db = 10*log10(Px_ave);

if(sigtype==3)
   for nn = 1:TotN;
        den = 1-ar_coef(1)*exp(-j*2*pi*(nn-1)/(2*TotN))-ar_coef(2)*exp(-j*4*pi*(nn-1)/(2*TotN))-ar_coef(3)*exp(-j*6*pi*(nn-1)/(2*TotN));
        if(den == 0)
            den = 0.000001;
        end;
        psd_true(nn) = vari/(abs(den)^2);       
   end;    
   subplot(2,2,3);
   if(vari~=0)
      psd_db = 10*log10(psd_true); 
      plot(index_f,psd_db);
   end;
   xlabel('Frequency (units of \pi)');
   ylabel('Magnitude (dB)');
   title('True power spectrum','fontsize',13);
   subplot(2,2,4);  
else
   subplot('Position',[0.13 0.11 0.77 0.344]);
end;
plot(index_f,Px_ave_db(1:512));
xlabel('Frequency (units of \pi)');
ylabel('Magnitude (dB)');
title('The ensemble average','fontsize',13);