welch.m

考虑L的三个不同值：L＝256（3个数据段）

%welch 法
vari=1;
N=1000;
nsect=10;
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 nn = 1:N;
        sigdata(nn) = sin(nn*pi/4+seta1)+sin(nn*pi/5+seta2)+0.2*sin(nn*pi/3);
    end
over=0.03;
win=4;
L=64;
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);
sigtype=2;
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);