Matlab第四版

clear;
close all;
fid=fopen('external-sinewave.dat','r')
a=fread(fid,'int32')
a=a'
N=length(a);
dt=0.0005
Fs=2000;
t=0:dt:(N-1)*dt
subplot(3,1,1);
plot(t,a)
xlabel('时间/s');
ylabel('采样值');
title('原始信号(时域)');
grid on;
b=xcorr(a,'unbiased')
z=fft(a,N);
Pyy =z.*conj(z)*2/N;
pyy=log(Pyy)
f =Fs*(0:(N-1))/N;
subplot(3,1,2)
plot(f,pyy(1:4096),'r')
xlabel('频率 (Hz)')
ylabel('能量')
title('频谱图(频域)')
grid on;
[c d]=max(Pyy)
f1=f(d)
k=0;
q=floor((N/2)/d);
k=1;
j=1;
for i=1:q
    Fpyy(k)=Pyy(i*d-2);
    Fpyy(k+1)=Pyy(i*d-1);
    Fpyy(k+2)=Pyy(i*d);
    Fpyy(k+3)=Pyy(i*d+1);
    Fpyy(k+4)=Pyy(i*d+2);
    z=[Fpyy(k) Fpyy(k+1) Fpyy(k+2) Fpyy(k+3) Fpyy(k+4)]
XB(j)=max(z);
k=k+5;
j=j+1;
end
thd=sqrt((sum(XB)-XB(1))/XB(1))
range='half';
noverlap=20;
window=hann(N);
[Pxx,f]=pwelch(a,window,noverlap,N,Fs,range);
Pxx1=10*log10(Pxx)
subplot(3,1,3)
plot(f,Pxx1)
fid2=fopen('result.txt','w');
fprintf(fid2,'\n数据基本信息：\r\n')
fprintf(fid2,'        采样点数 = %7.0f \r\n',length(a))      
fprintf(fid2,'        采样时间 = %7.3f s\r\n',length(a)*dt)                   
fprintf(fid2,'        采样频率 = %7.1f Hz\r\n',Fs)
fprintf(fid2,'        标准方差 = %7.3f \r\n',std(a))
fprintf(fid2,'         协 方 差 = %7.3f \r\n',cov(a))
fprintf(fid2,'     自相关系数 = %7.3f \r\n',corrcoef(a)) 

fprintf(fid2,'\n傅立叶变换结果：\r\n')
fprintf(fid2,'           基波频率= %1.3f Hz\r\n',f1)            
fprintf(fid2,'           基波周期 = %1.3f s\r\n',1/f(d)) 
fprintf(fid2,'           谐波畸变 = %1.3e \r\n',thd)