welch.m

there is a simple demo for non parameteric spectral estimation methods

function Px = welch(x)

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

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

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

HH = findobj(gcf,'Tag','wel_ol');
p = str2num(get(HH,'String'));
if (p >= 1) | (p < 0)  
   error('Overlap is invalid');
else
   over = p;
end;

HH = findobj(gcf,'Tag','wel_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','Welch 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 sigcnt=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((sigcnt-1)*N+nn);
    end;
elseif(sigtype == 3)
    sigdata(1) = noise((sigcnt-1)*N+1);
    sigdata(2) = ar_coef(1)*sigdata(1)+noise((sigcnt-1)*N+2);
    sigdata(3) = ar_coef(1)*sigdata(2)+ar_coef(2)*sigdata(1)+noise((sigcnt-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((sigcnt-1)*N+nn);
    end;
else
    for nn = 1:N;
        sigdata(nn) = noise((sigcnt-1)*N+nn);
    end;
end;

n_1 = 1;
n_0 = (1-over)*L;
nsect=1+floor((N-L)/(n_0));
Px=0;
for i=1:nsect
    Px_mper = mper(sigdata,win,n_1,n_1+L-1);
    Px = Px + Px_mper/nsect;
    n_1 = floor(n_1 + n_0);
end;
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 Welch 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);