%随机信号谱分析
clear
clc
close all hidden
%fni=input('随机信号谱分析－输入数据文件名','s');
%fid=fopen(fni,'r')
%分析内容（1=自谱,2=互谱,3=频响函数,4=相干函数）
%fun=fscanf(fid,'%f',1);
%sf=fscanf(fid,'%f',1); %读入采样频率值
%nfft=fscanf(fid,'%d',1); %读入FFT长度
%窗函数（1=矩形,2=汉宁,3=海明,4=布莱克曼,5=三角）
%win=fscanf(fid,'%d',1);
%fno=fscanf(fid,'%d',1); %读入输出数据文件名
%按行读入数据，第一行激励，第二行响应
%a=fscanf(fid,'%f',[2 inf]);
%status=fclose(fid);
fun=3;sf=500;nfft=2048;win=1;fno='out6_1.mat';
load y
a(1,:)=y;a(2,:)=y.*sin(y);
x=a(1,:);
y=a(2,:)-a(1,:);
f=0:sf/nfft:sf/2-sf/nfft;
switch win
case 1
w=boxcar(nfft);
case 2
w=hanning(nfft);
case 3
w=hamming(nfft);
case 4
w=blackman(nfft);
case 5
w=triang(nfft);
otherwise
w=boxcar(nfft);
end
switch fun
case 1
z=psd(y,nfft,sf,w,nfft/2);
case 2
z=csd(x,y,nfft,sf,w,nfft/2);
case 3
z=tfe(x,y,nfft,sf,w,nfft/2);
case 4
z=cohere(x,y,nfft,sf,w,nfft/2);
otherwise
;
end
%幅频曲线
nn=1:nfft/4;
subplot(2,1,1);
plot(f(nn),abs(z(nn)));
xlabel('频率(Hz)');
ylabel('幅值');
grid on;
if fun>1&fun<4
%相频曲线
subplot(2,1,2);
plot(f(nn),angle(z(nn)));
xlabel('频率(Hz)');
ylabel('相位');
grid on;
end
fid=fopen(fno,'w');
if fun>1&fun<4
for k=1:nfft/2
fprintf(fid,'%f%f\n',f(k),abs(z(k)));%自谱和相干数据
end
else
for k=1:nfft/2
fprintf(fid,'%f%f%f\n',f(k),real(z(k)),imag(z(k)));%互谱和频响函数数据
end
end
status=fclose(fid);
