function moments = spectral_moments(signal, fs, type);
% siglen=length(signal);
siglen=2048;
switch type
    case 'fft'
        spect=fft(signal,2048); %twosided spectrum .*hamming(sig_len)
        spect=abs(spect);
        spect=spect(1:floor(siglen/2));
        spect(2:end)=2*spect(2:end);
        spect=spect.^2/siglen; %power spectrum
        f = (fs*(0:(siglen/2)-1)/siglen)';
        norm_spect=spect/(sum(spect));
        m1=sum(f.*norm_spect);
        m2=sqrt(sum(((f-m1).^2).*norm_spect));
        m3=(sum(((f-m1).^3).*norm_spect))/m2^3;
        m4=(sum(((f-m1).^4).*norm_spect))/m2^4;  
       
        %slope
        p=polyfit(f(1:find(f>2000,1)),spect(1:find(f>2000,1)),1);
        sl=p(1);
        moments=[m1 m2 m3 m4 sl];
    case 'welch'
        [Pxx,F] = pwelch(signal,[],[],siglen-1,fs);
%         HdB=10*log10(Pxx); Pxx=HdB-max(HdB);
        norm_spect=Pxx/(sum(Pxx));
        m1=sum(F.*norm_spect);
        m2=sqrt(sum(((F-m1).^2).*norm_spect));
        m3=(sum(((F-m1).^3).*norm_spect))/m2^3;
        m4=(sum(((F-m1).^4).*norm_spect))/m2^4;  
        p=polyfit(F(1:find(F>2000,1)),spect(1:find(F>2000,1)),1);
        sl=p(1);
        moments=[m1 m2 m3 m4 sl];
    case 'multitaper'
        [pxx,f] = pmtm(signal,3,length(signal),fs);
        norm_spect=pxx/sum(pxx);
        m1=sum(f.*norm_spect);
        m2=sqrt(sum(((f-m1).^2).*norm_spect));
        m3=(sum(((f-m1).^3).*norm_spect))/m2^3;
        m4=(sum(((f-m1).^4).*norm_spect))/m2^4;  
        p=polyfit(f(1:find(f>2000,1)),spect(1:find(f>2000,1)),1);
        sl=p(1);
        moments=[m1 m2 m3 m4 sl];
    case 'lpc'
        p=20;
        [a, g]= lpc(signal,p);
        [H, f]=freqz(sqrt(g),a,1024,fs);
        lpcsp=((abs(H).^2));
%         lpcsp=10*log10((abs(H).^2));
        norm_spect=lpcsp/sum(lpcsp);
        m1=sum(f.*norm_spect);
        m2=sqrt(sum(((f-m1).^2).*norm_spect));
        m3=(sum(((f-m1).^3).*norm_spect))/m2^3;
        m4=(sum(((f-m1).^4).*norm_spect))/m2^4;  
        p=polyfit(f(1:find(f>2000,1)),spect(1:find(f>2000,1)),1);
        sl=p(1);
        moments=[m1 m2 m3 m4 sl];
end