glp.m

hosa包

clear all;
%% open data file 
fid=fopen('d:\HKWH_Z1.ABF','r');%file path
xdata=fread(fid,'float32');   %read data from file
fclose(fid);%close file
clear fid;

%% extract 32768 points from the data file 
i=1;  %% change to the own number for each student i=9
startp=i*120000; %% start point of the analyzed data
endp=startp+32767; %% end point of the analyzed data
x1=xdata(startp:endp); %% data extraction

fs=4000; %% sampling frequency in Hz
ts=1/fs; %% sampling periodic in second

tx=0:ts:(1023*ts); %% vector for time achse

figure(1); 
%% plot the long time signal
subplot(211);
plot(xdata);
xlabel('sampling point');
ylabel('vibration amplitude/um');

%% plot the extracted signal
subplot(212);
plot(tx,x1(1:1024));
xlabel('time / second');
ylabel('vibration amplitude / um');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%based on Welch's method power spectral density estimate
[Pxx,F]=pwelch(x1,[],[],[],fs,'onesided');% the pwelch function
figure(2);
%plot power spectral density
semilogy(F(1:300),Pxx(1:300));
xlabel('frequency / Hz');%x label
ylabel('APSD /  um**2/Hz');%y label
title('Power Spectral Density')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%双谱估计
figure(3);
[bspec,waxis]=bispecd(x1,16);
title('双谱估计平面图')

figure(4);
x2=abs(bspec);
mesh(x2)
title('双谱估计三维图')