multitap.m

工程计算MATLAB code to calculate the reorthogonalized sine tapers input: N = the length of the time se

function [f, s] = multitap(x, taper)

%%[f, s]=multitap(x, taper)
%%x a data vector, of length N; taper a matrix (N-by-p) of p tapers
%%each of length N, as the data to be tapered
%%returns the frequency f normalized to be in [0, 0.5],  and the
%%multitapered spectrum estimate s, as the average of the p individual
%%eigenspectra.  the individual eigenspectra are just the modulus
%%square of the fft of the data x multiplied by one of the tapers.
%%plots the spectrum on Decibel scale (10*log10(s)) vs. frequency
%%no taper at all, the periodogram, corresponds to taper=ones(N,1);


t=size(taper);
N=t(1);
p=t(2);

s=zeros(N/2,1 );

for i=1:p
	y=x.*taper(:,i);
	[freq, a]=spec(y, N);
	s = s+a;
end
s=s/p;
f=freq;
% figure
% plot(f, 10*log10(s)), title('Multitapered Spectrum Estimate')
% xlabel('Frequency (cycles/unit time)'), ylabel('Power (dB)')