另外,还有为啥得到的谱图只有到25Hz?
采用频率为100.16Hz,FFT长度1024
程序如下:
%随机信号谱分析
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%
fni=input('随机信号谱分析-输入数据文件名:','s');
fid=fopen(fni,'r');
fun = fscanf(fid,'%d',1); %分析内容(1自谱,2互谱,3频响函数,4相干函数)
sf = fscanf(fid,'%f',1); %采样频率
nfft = fscanf(fid,'%d',1); %FFT长度
win = fscanf(fid,'%d',1); %窗函数(1矩形,2汉宁,3海明,4布莱克曼,5三角)
fg = fscanf(fid,'%d',1); %图形方式(1实部、虚部图,2幅值、相位图)
fno = fscanf(fid,'%s',1); %输出数据文件名
a = fscanf(fid,'%f',[2,inf]); %按列读入数据,第1列激励,第2列响应
x=a(1,:);
y=a(2,:)-a(1,:);
status=fclose(fid);
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);
if fg==1
plot(f(nn),real(z(nn)));
ylabel('实部'); %添加纵向坐标轴的标注
else
plot(f(nn),abs(z(nn)));
ylabel('幅值');
end
xlabel('频率 (Hz)'); %添加横向坐标轴的标注
grid on; %在图幅上添加坐标网格
if fun>1&fun<4
subplot(2,1,2);
if fg==1
plot(f(nn),imag(z(nn)));
ylabel('虚部'); %添加纵向坐标轴的标注
else
plot(f(nn),angle(z(nn)));
ylabel('相位');
end
xlabel('频率 (Hz)'); %添加横向坐标轴的标注
grid on; %在图幅上添加坐标网格
end
fid=fopen(fno,'w'); %以写的方式打开文件或建立一个新文件
if fun==1|fun==4 %自谱和相干函数
for k=1:nfft/2
fprintf(fid,'%f %f\r\n',f(k),abs(z(k))); %每行输出2个实型数据,f为频率,z为幅值
end
else %互谱和频响函数
for k=1:nfft/2
fprintf(fid,'%f %f %f\r\n',f(k),real(z(k)),imag(z(k))); %每行输出3个实型数据
end
end
status=fclose(fid);
|