Matlab第三版

clear;
close all;
fid=fopen('external-sinewave.dat','r')
a=fread(fid,'int32')
a=a'
t=0:0.0005:4095*0.0005
 
subplot(3,1,1);
plot(t,a)
xlabel('时间/s');
ylabel('采样值');
title('原始信号(时域)');
grid on;
%a=a-mean(a);
Fs=2000;
N=4096;
z=fft(a);
Pyy =z.*conj(z)/N;
pyy=log(Pyy)
j=mean(pyy,2);
f = 2000*(0:4095)/N;
subplot(3,1,2)
plot(f,pyy(1:4096),'r')
k=0;
for i=2:2048
    if pyy(i)-pyy(i-1)>2 &&pyy(i)-pyy(i+1)>2
        k=k+1;
fpyy(k)=pyy(i);
Fpyy(k)=Pyy(i);
end
end

xlabel('频率 (Hz)')
ylabel('能量')
title('频谱图(频域)')
grid on;

[b c]=max(Pyy(1:N));
o=2048/c
p=floor(o)
g=0;
for q=1:p
    g=g+1;
    ffpyy(g)=Pyy(q*c-1);
end
fffpyy=log(ffpyy)
h = 15.625*(0:53);

subplot(3,1,3)
plot(h,fpyy(1:54),'r') 
xlabel('频率 (Hz)')
ylabel('能量')
title('频谱图(频域)')
grid on;
sum=sum(Fpyy)
maxf=max(Fpyy)
thd=(sum-maxf)/maxf

   
fprintf('\n数据基本信息：\n')
fprintf('        采样点数 = %7.0f \n',length(a))      
fprintf('        采样时间 = %7.3f ms\n',length(a)*2)                   
fprintf('        采样频率 = %7.1f Hz\n',1000/0.5)
fprintf('        标准方差 = %7.3f \n',std(a))
fprintf('         协 方 差 = %7.3f \n',cov(a))
fprintf('     自相关系数 = %7.3f \n\n',corrcoef(a)) 

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