%参考文献上采用的multitaper的方法，
Fs=2000;
nw=3; %时间带宽积
N = length(x);
%nfft=2048；
%[Pxx,f]= pmtm(signal,nw,nfft,samplerate);
[Pxx,f]= pmtm(signal,nw,[0:250],Fs);%这种是默认的nfft
plot_Pxx=10*log10(Pxx);
figure(2);
plot(f,plot_Pxx);
set(gca,'xlim',[0 250]);
%NMPK工具包里和neuroexplorer说明书里的采用Periodogram Using FFT方法
%NMPK工具包里是NSxPowerSpectrum(NS3, channelNumber)函数
x=double(NS3.Data(1,:));
samplerate = 2000;
N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(samplerate*N)).*abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:samplerate/length(x):samplerate/2;
plot(freq,10*log10(psdx), 'b'); grid on;
title('Periodogram Using FFT');
xlabel('Frequency (Hz)'); ylabel('Power/Frequency (dB/Hz)');
set(gca,'xlim',[0 200]);
%welch窗函数法
x=double(NS3.Data(1,:));
samplerate = 2000;
N = length(x);
%nfft=max(2^(nextpow2(N)+0),N);;%%
window=hamming(100);
noverlap=20;
range='onesided';
%[Pxx1,f]=pwelch(x,window,noverlap,N,samplerate,range);
[Pxx1,f]=pwelch(x,window,noverlap,2048,samplerate,range)
plot_Pxx1=10*log10(Pxx1);
plot(f,plot_Pxx1);
title('Welch Hamming');
ylabel('db/hz');
%set(gca,'xlim',[0 200]);
chronux工具包mtspectrumc函数
openNSx;
data=double(NS3.Data(1,:));%只提取第一行
%data=double(NS3.Data);%提取所有行
matlab_params;%调出param脚本
movingwin=[0.5 0.05]; % set the moving window dimensions
%data= ;
[S,f]=mtspectrumc(data,params);
plot_vector(S,f,[],[],'m');
set(gca,'FontName','Times New Roman','Fontsize', 16);
legend(['W=' num2str(3/size(data,1)*params.samplerate) 'Hz'],['W=' num2str(params.tapers(1)/size(data,1)*params.samplerate) 'Hz']);
title('LFP spectrum for 1 trial with varying bandwidths');
position_figure
xlabel('Frequency/Hz');
评论0