admin 管理员组文章数量: 1087652
2024年3月7日发(作者：eof错误代码)
程序4-1
%最小二乘法消除多项式趋势项
%%%%%%%%%%%%%%%%%%%%%%%%
clear % 清除内存中所有变量和函数
clc % 清除工作窗口中所显示的内容
close all hidden % 关闭所有隐藏的窗口
%%%%%%%%%%%%%%%%%%%%%%%%
%提示用键盘输入输入数据文件名
fni=input('消除多项式趋势项-输入数据文件名:','s');
%以只读方式打开数据文件
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %读入采样频率值
m = fscanf(fid,'%d',1); %读入拟合多项式阶数
fno = fscanf(fid,'%s',1);%读入输出数据文件名
x = fscanf(fid,'%f',inf);%读入时程数据存成列向量
%关闭数据文件
status=fclose(fid);
%取信号数据长度
n=length(x);
%建立离散时间列向量
t=(0:1/sf:(n-1)/sf)';
%计算趋势项的多项式待定系数向量a
a=polyfit(t,x,m);
%用x减去多项式系数a生成的趋势项
y=x-polyval(a,t);
%将分成2行1列的图形窗口的第1列设为当前绘图区域
subplot(2,1,1);
%绘制x对于t的时程曲线图形
plot(t,x);
%在图幅上添加坐标网格
grid on;
%将分成2行1列的图形窗口的第2列设为当前绘图区域
subplot(2,1,2);
%绘制y对于t的时程曲线图形
plot(t,y);
%在图幅上添加坐标网格
grid on;
%以写的方式打开文件或建立一个新文件
fid=fopen(fno,'w');
%进行n次循环将计算结果写到输出数据文件中
for k=1:n
%每行输出两个实型数据，t为时间，y为消除趋势项后的结果
fprintf(fid,'%f %fn',t(k),y(k));
%循环体结束语句
end
%关闭数据文件
status=fclose(fid);
程序4-2
%五点滑动平均法平滑处理
%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%
fni=input('五点滑动平均法平滑处理-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
m = fscanf(fid,'%d',1); %平滑次数
fno = fscanf(fid,'%s',1);%输出数据文件名
x = fscanf(fid,'%f',inf);%输入数据存成列向量
status=fclose(fid);
%取信号数据长度
n=length(x);
%建立离散时间列向量
t=(0:1/sf:(n-1)/sf)';
%将x赋值给a
a=x;
%循环m次进行平滑处理计算
for k=1:m
b(1)=(3*a(1)+2*a(2)+a(3)-a(4))/5;
b(2)=(4*a(1)+3*a(2)+2*a(3)+a(4))/10;
for j=3:n-2
b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2))/5;
end
b(n-1)=(a(n-3)+2*a(n-2)+3*a(n-1)+4*a(n))/10;
b(n)=(-a(n-3)+a(n-2)+2*a(n-1)+3*a(n))/5;
a=b;
end
%将a赋值给y
y=a;
%将分成2行1列的图形窗口的第1列设为当前绘图区域
subplot(2,1,1);
%绘制平滑前的时程曲线图形
plot(t,x);
%添加横向坐标轴的标注
xlabel('时间 (s)');
%添加纵向坐标轴的标注
ylabel('加速度(g)');
%在图幅上添加坐标网格
grid on;
%将分成2行1列的图形窗口的第2列设为当前绘图区域
subplot(2,1,2);
%绘制平滑后的时程曲线图形
plot(t,y);
%添加横向坐标轴的标注
xlabel('时间 (s)');
%添加纵向坐标轴的标注
ylabel('加速度(g)');
%在图幅上添加坐标网格
grid on;
%打开文件输出平滑后的数据
fid=fopen(fno,'w');
for k=1:n
%每行写两个实型数据，t为时间，y为平滑后的数据
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序4-2
%五点滑动平均法平滑处理
%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%
fni=input('五点滑动平均法平滑处理-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
m = fscanf(fid,'%d',1); %平滑次数
fno = fscanf(fid,'%s',1);%输出数据文件名
x = fscanf(fid,'%f',inf);%输入数据存成列向量
status=fclose(fid);
%取信号数据长度
n=length(x);
%建立离散时间列向量
t=(0:1/sf:(n-1)/sf)';
%将x赋值给a
a=x;
%循环m次进行平滑处理计算
for k=1:m
b(1)=(3*a(1)+2*a(2)+a(3)-a(4))/5;
b(2)=(4*a(1)+3*a(2)+2*a(3)+a(4))/10;
for j=3:n-2
b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2))/5;
end
b(n-1)=(a(n-3)+2*a(n-2)+3*a(n-1)+4*a(n))/10;
b(n)=(-a(n-3)+a(n-2)+2*a(n-1)+3*a(n))/5;
a=b;
end
%将a赋值给y
y=a;
%将分成2行1列的图形窗口的第1列设为当前绘图区域
subplot(2,1,1);
%绘制平滑前的时程曲线图形
plot(t,x);
%添加横向坐标轴的标注
xlabel('时间 (s)');
%添加纵向坐标轴的标注
ylabel('加速度(g)');
%在图幅上添加坐标网格
grid on;
%将分成2行1列的图形窗口的第2列设为当前绘图区域
subplot(2,1,2);
%绘制平滑后的时程曲线图形
plot(t,y);
%添加横向坐标轴的标注
xlabel('时间 (s)');
%添加纵向坐标轴的标注
ylabel('加速度(g)');
%在图幅上添加坐标网格
grid on;
%打开文件输出平滑后的数据
fid=fopen(fno,'w');
for k=1:n
%每行写两个实型数据，t为时间，y为平滑后的数据
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序4-3
%滑动平均法消除趋势项
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%
fni=input('滑动平均法消除趋势项-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
l = fscanf(fid,'%d',1); %滑动阶次
m = fscanf(fid,'%d',1); %平滑次数
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
%取信号数据长度
n=length(x);
%建立离散时间列向量
t=(0:1/sf:(n-1)/sf);
%生成一个元素全为1的行向量
b=ones(1,l);
%信号两端分别向外延伸l个数据
a=[b*x(1),x,b*x(n)];
b=a;
%按平滑次数循环进行滑动平均处理计算趋势项
for k=1:m
for j=l+1:n-l
b(j)=mean(a(j-l:j+l));
end
a=b;
end
%用输入信号x减去与平滑趋势项a
y=x(1:n)-a(l+1:n+l);
%同时绘制x对于t和y对于t的时程曲线
plot(t,x,':',t,y,t,a(l+1:n+l),'-.');
%添加横向坐标轴的标注
xlabel('时间 (s)');
%添加纵向坐标轴的标注
ylabel('位移 mm');
%在图幅上添加图例
legend('输入','输出','趋势');
%在图幅上添加坐标网格
grid on;
%以写的方式打开文件或建立一个新文件
fid=fopen(fno,'w');
%进行n次循环将结果写到输出数据文件中
for k=1:n
%每行写两个实型数据，t为时间，y为消除趋势项后的结果
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序4-4
%五点三次法平滑处理(时域和频域)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('五点三次平滑处理-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
it = fscanf(fid,'%d',1); %数据类型(1时域,2频域)
m = fscanf(fid,'%d',1); %平滑次数
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[it,inf]);%输入数据存成行向量
status=fclose(fid);
%取信号数据长度
n=length(x(1,:));
for l=1:it
a=x(l,:);
%循环m次进行平滑处理
for k=1:m
b(1)=(69*a(1)+4*(a(2)+a(4))-6*a(3)-a(5))/70;
b(2)=(2*(a(1)+a(5))+27*a(2)+12*a(3)-8*a(4))/35;
for j=3:n-2
b(j)=(-3*(a(j-2)+a(j+2))+12*(a(j-1)+a(j+1))+17*a(j))/35;
end
b(n-1)=(2*(a(n)+a(n-4))+27*a(n-1)+12*a(n-2)-8*a(n-3))/35;
b(n)=(69*a(n)+4*(a(n-1)+a(n-3))-6*a(n-2)-a(n-4))/70;
a=b;
end
y(l,:)=a;
end
%绘制平滑前后的曲线图形
if it==1 %时域信号
%建立离散时间向量
nn=1:2000;
t=0:1/sf:(n-1)/sf;
%同时绘制x对于t和y对于t的时域曲线
plot(t(nn),x(nn),':',t(nn),y(nn));
xlabel('时间 (s)'); %添加横向坐标轴的标注
ylabel('幅值'); %添加纵向坐标轴的标注
legend('平滑前','平滑后');%在图幅上添加图例
grid on;
else %频域信号
%建立离散频率向量
nn=1:256;
f=0:sf/n:(n-1)*sf/n;
%同时绘制x的实部对于f和y的实部对于f的频域曲线
subplot(2,1,1);
plot(f(nn),x(1,nn),':',f(nn),y(1,nn));
xlabel('频率 (Hz)'); %添加横向坐标轴的标注
ylabel('实部'); %添加纵向坐标轴的标注
legend('平滑前','平滑后');%在图幅上添加图例
grid on;
%同时绘制x的虚部对于f和y的虚部对于f的频域曲线
subplot(2,1,2);
plot(f(nn),x(2,nn),':',f(nn),y(2,nn));
xlabel('频率 (Hz)'); %添加横向坐标轴的标注
ylabel('虚部'); %添加纵向坐标轴的标注
legend('平滑前','平滑后');%在图幅上添加图例
grid on;
end
%打开文件输出平滑后的数据
fid=fopen(fno,'w');
for k=1:n
if it==1 %时域信号输出时间和幅值
fprintf(fid,'%f %fn',t(k),y(k));
else %频域信号输出频率、实部和虚部
fprintf(fid,'%f %f %fn',f(k),y(1,k),y(2,k));
end
end
status=fclose(fid);
程序5-1
%频域低通和带通滤波
%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%
fni=input('频域带通滤波-输入数据文件名:','s');
fid = fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
fmin = fscanf(fid,'%f',1); %最小截止频率
fmax = fscanf(fid,'%f',1); %最大截止频率
sx = fscanf(fid,'%s',1); %读入横向坐标轴的标注
sy = fscanf(fid,'%s',1); %读入纵向坐标轴的标注
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
%取信号数据长度
n=length(x);
%建立离散时间列向量
t=(0:1/sf:(n-1)/sf)';
%取大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%四舍五入取整求最小截止频率对应数组元素的下标
ni=round(fmin*nfft/sf+1);
%四舍五入取整求最大截止频率对应数组元素的下标
na=round(fmax*nfft/sf+1);
%进行FFT变换,结果存于y
y=fft(x,nfft);
%建立一个长度为nfft元素全为0的向量
a=zeros(1,nfft);
%将y的正频率带通内的元素赋值给a
a(ni:na)=y(ni:na);
%将y的负频率带通内的元素赋值给a
a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
%进行FFT逆变换,结果存于y
y=ifft(a,nfft);
%取逆变换的实部n个元素为滤波结果列向量
y=(real(y(1:n)))';
%绘制滤波前的时程曲线图形
subplot(2,1,1);
plot(t,x);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy);
grid on;
%绘制滤波后的时程曲线图形
subplot(2,1,2);
plot(t,y);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy);
grid on;
%打开文件输出滤波后的数据
fid=fopen(fno,'w');
for k=1:n
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序5-2
%频域高通和带阻滤波
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%
fni=input('频域带阻滤波-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
fmin = fscanf(fid,'%f',1); %最小截止频率
fmax = fscanf(fid,'%f',1); %最大截止频率
sx = fscanf(fid,'%s',1); %读入横向坐标轴的标注
sy = fscanf(fid,'%s',1); %读入纵向坐标轴的标注
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
%取信号数据长度
n=length(x);
%建立离散时间列向量
t=(0:1/sf:(n-1)/sf)';
%取大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%四舍五入取整求最小截止频率对应数组元素的下标
ni=round(fmin*nfft/sf+1);
%四舍五入取整求最大截止频率对应数组元素的下标
na=round(fmax*nfft/sf+1);
%进行FFT变换,结果存于y
y=fft(x,nfft);
%建立一个长度为nfft元素全为0的向量
a=zeros(1,nfft);
%将y的正频率带阻内的元素赋值为0
y(ni:na)=a(ni:na);
%将y的负频率带阻内的元素赋值为0
y(nfft-na+1:nfft-ni+1)=a(nfft-na+1:nfft-ni+1);
%进行IFFT变换,结果存于a
a=ifft(y,nfft);
%取逆变换的实部n个元素为滤波结果列向量
y=real(a(1:n));
%绘制滤波前的时程曲线图形
subplot(2,1,1);
plot(t,x);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy);
grid on;
%绘制滤波后的时程曲线图形
subplot(2,1,2);
plot(t,y);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy);
grid on;
%打开文件输出滤波后的数据
fid=fopen(fno,'w');
for k=1:n
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序5-3
%运用完全设计法的IIR滤波器滤波
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('运用完全设计法的IIR滤波器滤波-输入数据文件名:','s');
fid=fopen(fni,'r');
fs = fscanf(fid,'%f',1); %采样频率
%滤波方式选择（1低通,2高通,3带通,4带阻）
fun = fscanf(fid,'%d',1);
%滤波器选择（1巴特沃斯,2切比雪夫Ⅰ,3切比雪夫Ⅱ,4椭圆）
mod = fscanf(fid,'%d',1);
if fun<=2
wp = fscanf(fid,'%f',1);%通带截止频率(Hz)
ws = fscanf(fid,'%f',1);%阻带截止频率(Hz)
else
wp = fscanf(fid,'%f',2);%通带截止频率(Hz)
ws = fscanf(fid,'%f',2);%阻带截止频率(Hz)
end
rp = fscanf(fid,'%f',1); %通带波动系数(dB)
rs = fscanf(fid,'%f',1); %阻带衰减系数(dB)
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
%确定滤波方式
switch fun
case 1 %低通
ft='low';
case 2 %高通
ft='high';
case 3 %带通
ft='bandpass';
case 4 %带阻
ft='stop';
otherwise
ft='low';
end
%根据滤波器种类进行IIR滤波器设计
switch mod
%巴特沃斯滤波器
case 1
[n wn]=buttord(wp/(fs/2),ws/(fs/2),rp,rs);
[b a]= butter(n,wn,ft);
%切比雪夫Ⅰ型滤波器
case 2
[n wn]=cheb1ord(wp/(fs/2),ws/(fs/2),rp,rs);
[b a]= cheby1(n,rp,wn,ft);
%切比雪夫Ⅱ型滤波器
case 3
[n wn]=cheb2ord(wp/(fs/2),ws/(fs/2),rp,rs);
[b a]= cheby2(n,rs,wn,ft);
%椭圆滤波器
case 4
[n wn]=ellipbuttord(wp/(fs/2),ws/(fs/2),rp,rs);
[b a]= ellip(n,rp,rs,wn,ft);
end
%计算滤波器的频率响应
[h w]=freqz(b,a,1024,fs);
m=length(x);
t=0:1/fs:(m-1)/fs;
%用设计出的滤波器进行滤波
y=filter(b,a,x);
%绘制滤波器的频率响应图
subplot(2,1,1);
plot(w,abs(h));
%添加横向坐标轴的标注
xlabel('频率 (Hz)');
%添加纵向坐标轴的标注
ylabel('幅值');
%在图幅上添加坐标网格
grid on;
%绘制滤波前后的信号时程图
subplot(2,1,2);
plot(t,x,':',t,y);
%添加横向坐标轴的标注
xlabel('时间 (s)');
%添加纵向坐标轴的标注
ylabel('加速度(g)');
%在图幅上添加图例
legend('滤波前','滤波后');
grid on;
%打开文件输出滤波后的数据
fid=fopen(fno,'w');
for k=1:m
fprintf(fid,'%f %f n',t(k),y(k));
end
status=fclose(fid);
程序5-4
%多带FIR滤波器滤波
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('多带FIR滤波器滤波-输入数据文件名:','s');
fid=fopen(fni,'r');
fs = fscanf(fid,'%f',1); %采样频率
%滤波器函数选择（1=fir2,2=firls,3=remez,4=cremez）
mod = fscanf(fid,'%d',1);
n = fscanf(fid,'%d',1); %滤波器阶数
k = fscanf(fid,'%d',1); %频率点数
f = fscanf(fid,'%f',k); %指定频率点(Hz)
%通阻状态(0=阻 or 1=通)，对于函数cremez为期望幅值
a = fscanf(fid,'%f',k);
w = fscanf(fid,'%f',k/2); %权向量
fno = fscanf(fid,'%s',1); %输出数据文件名
%输入数据存成行向量（1加噪信号，2原始信号）
x = fscanf(fid,'%f',[2,inf]);
status=fclose(fid);
m=length(x(1,:));
t=0:1/fs:(m-1)/fs;
%频率归一化处理
f=2*f/fs;
%多带FIR滤波器设计
switch mod
case 1 %fir2滤波函数
b= fir2(n,f,a);
case 2 %firls滤波函数
b= firls(n,f,a,w);
case 3 %remez滤波函数
b= remez(n,f,a,w);
case 4 %cremez滤波函数
b= cremez(n,f,a,w,'real');
end
%用多带FIR滤波器进行滤波
v=filter(b,1,x(1,:));
%消除延时处理
y=zeros(1,m);
y(1:m-n/2)=v(n/2+1:m);
%计算滤波器的频率响应
[h w]=freqz(b,1,1024,fs);
%绘制滤波器的频率响应图
subplot(2,1,1)
plot(w,abs(h));
%添加横向坐标轴的标注
xlabel('频率 (Hz)');
%添加纵向坐标轴的标注
ylabel('幅值');
grid on;
%绘制滤波前后的信号时程比较图
subplot(2,1,2)
nn=1:200;
plot(t(nn),x(1,nn),'-.',t(nn),x(2,nn),':',t(nn),y(nn));
%添加横向坐标轴的标注
xlabel('时间 (s)');
%添加纵向坐标轴的标注
ylabel('幅值');
%在图幅上添加图例
legend('加噪信号','原始信号','滤波信号');
grid on;
%打开文件输出滤波后的数据
fid=fopen(fno,'w');
for k=1:m
fprintf(fid,'%f %f n',t(k),y(k));
end
status=fclose(fid);
程序5-5
%频域积分
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%
fni=input('频域积分-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
fmin = fscanf(fid,'%f',1); %最小截止频率
fmax = fscanf(fid,'%f',1); %最大截止频率
c = fscanf(fid,'%f',1); %单位变换系数
it = fscanf(fid,'%f',1); %积分次数
sx = fscanf(fid,'%s',1); %横向坐标轴的标注
sy1 = fscanf(fid,'%s',1); %纵向坐标轴输入单位的标注
sy2 = fscanf(fid,'%s',1); %纵向坐标轴输出单位的标注
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
n=length(x);
%建立时间向量
t=0:1/sf:(n-1)/sf;
%大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%FFT变换
y=fft(x,nfft);
%计算频率间隔（Hz/s）
df=sf/nfft;
%计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);
%计算圆频率间隔（rad/s）
dw=2*pi*df;
%建立正的离散圆频率向量
w1=0:dw:2*pi*(0.5*sf-df);
%建立负的离散圆频率向量
w2=2*pi*(0.5*sf-df):-dw:0;
%将正负圆频率向量组合成一个向量
w=[w1,w2];
%以积分次数为指数，建立圆频率变量向量
w=w.^it;
%进行积分的频域变换
a=zeros(1,nfft);
a(2:nfft-1)=y(2:nfft-1)./w(2:nfft-1);
if it==2
%进行二次积分的相位变换
y=-a;
else
%进行一次积分的相位变换
real(y)=imag(a);
imag(y)=-real(a);
end
a=zeros(1,nfft);
%消除指定正频带外的频率成分
a(ni:na)=y(ni:na);
%消除指定负频带外的频率成分
a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
%IFFT变换
y=ifft(a,nfft);
%取逆变换的实部n个元素并乘以单位变换系数为积分结果
y=real(y(1:n))*c;
%绘制积分前的时程曲线图形
subplot(2,1,1);
plot(t,x);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy1);
grid on;
%绘制积分后的时程曲线图形
subplot(2,1,2);
plot(t,y);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy2);
grid on;
%打开文件输出积分后的数据
fid=fopen(fno,'w');
for k=1:n
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序5-6
%频域微分
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%
fni=input('频域微分-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
fmin = fscanf(fid,'%f',1); %最小截止频率
fmax = fscanf(fid,'%f',1); %最大截止频率
c = fscanf(fid,'%f',1); %单位变换系数
it = fscanf(fid,'%f',1); %微分次数
sx = fscanf(fid,'%s',1); %横向坐标轴的标注
sy1 = fscanf(fid,'%s',1); %纵向坐标轴输入单位的标注
sy2 = fscanf(fid,'%s',1); %纵向坐标轴输出单位的标注
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
n=length(x);
%建立时间向量
t=0:1/sf:(n-1)/sf;
%大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%FFT变换
y=fft(x,nfft);
%计算频率间隔（Hz/s）
df=sf/nfft;
%计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);
%计算圆频率间隔（rad/s）
dw=2*pi*df;
%建立从负到正的圆频率向量
w1=0:dw:2*pi*(0.5*sf-df);
w2=2*pi*(0.5*sf-df):-dw:0;
w=[w1,w2];
%以微分次数为指数，建立圆频率变量向量
w=w.^it;
%进行微分的频域变换
a=y.*w;
if it==2
%进行二次微分的相位变换
y=-a;
else
%进行一次微分的相位变换
real(y)=imag(a);
imag(y)=-real(a);
end
a=zeros(1,nfft);
%消除指定正频带外的频率成分
a(ni:na)=y(ni:na);
%消除指定负频带外的频率成分
a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
%IFFT变换
y=ifft(a,nfft);
%取逆变换的实部n个元素并乘以单位变换系数为微分结果
y=real(y(1:n))*c;
%绘制微分前的时程曲线图形
subplot(2,1,1);
plot(t,x);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy1);
grid on;
%绘制微分后的时程曲线图形
subplot(2,1,2);
plot(t,y);
%添加横向坐标轴的标注
xlabel(sx);
%添加纵向坐标轴的标注
ylabel(sy2);
grid on;
%打开文件输出微分后的数据
fid=fopen(fno,'w');
for k=1:n
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序5-7
%相关函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%
fni=input('相关函数-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
%分析内容(1=自相关,2=互相关)
fun = fscanf(fid,'%d',1);
m = fscanf(fid,'%d',1); %相关数据长度（点数）
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[fun,inf]);%输入数据存成行向量
status=fclose(fid);
n=length(x(1,:));
%建立离散时间向量
t=0:1/sf:(n-1)/sf;
if fun==1
%计算自相关函数向量
a=xcorr(x(1,:),m);
else
%计算互相关函数向量
a=xcorr(x(1,:),x(2,:),m);
end
%取正频率段的相关函数
y=a(m+1:2*m+1);
%绘制输入数据的时程曲线图形
subplot(2,1,1);
nn=1:4000;
if fun==1
plot(t(nn),x(1,nn));
else
plot(t(nn),x(1,nn),':',t(nn),x(2,nn));
legend('x','y');
end
grid on;
%绘制相关函数曲线图形
subplot(2,1,2);
nn=1:m;
plot(t(nn),y(nn));
grid on;
%打开文件输出相关函数数据
fid=fopen(fno,'w');
for k=1:m
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序6-1
%随机信号谱分析
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%
fni=input('随机信号谱分析-输入数据文件名:','s');
fid=fopen(fni,'r');
%分析内容(1=自谱,2=互谱,3=频响函数,4=相干函数)
fun = fscanf(fid,'%d',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,'%s',1); %输出数据文件名
%按行读入数据,第1行激励,第2行响应
a = fscanf(fid,'%f',[2,inf]);
status=fclose(fid);
%取激励信号向量
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
%每行输出2个实型数据，f为频率，z为幅值
fprintf(fid,'%f %fn',f(k),abs(z(k)));
end
%输出互谱和频响函数数据
else
for k=1:nfft/2
%每行输出3个实型数据，f为频率，z的实、虚部
fprintf(fid,'%f %f %fn',f(k),real(z(k)),imag(z(k)));
end
end
status=fclose(fid);
程序6-2
%窗函数（MATLAB中没有这些窗函数）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('窗函数-输入数据文件名:','s');
fid=fopen(fni,'r');
n = fscanf(fid,'%d',1); %数据长度
%窗函数(1=余弦坡度窗,2=帕曾窗,3=指数窗,4=高斯窗)
win = fscanf(fid,'%d',1);
%参数(除指数类窗需要输入参数外其余输入1)
c = fscanf(fid,'%f',1);
fno = fscanf(fid,'%s',1);%输出数据文件名
status=fclose(fid);
%建立数据长度向量
l=0:n-1;
%计算窗函数
switch win
case 1 %余弦坡度窗
m=round(n*4/5);
w=[ones(1,m),(1+cos(5*pi*l(m:n-1)/n))/2];
case 2 %帕曾窗
w=[1-6*(l(1:n/2)/n).^2+6*(l(1:n/2)/n).^3,2*(1-l(n/2+1:n)/n).^3];
case 3 %指数窗
w=exp(-c*l);
case 4 %高斯窗
w=exp(-c*l.^2);
end
%绘制窗函数曲线图形
plot(l/n,w);
grid on;
%打开文件输出窗函数的数据
fid=fopen(fno,'w');
for k=1:n
fprintf(fid,'%d %fn',k,w(k));
end
status=fclose(fid);
程序6-3
%细化FFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%
fid=fopen(fni,'r');
fs = fscanf(fid,'%f',1); %采样频率
fi = fscanf(fid,'%f',1); %最小细化截止频率
np = fscanf(fid,'%d',1); %放大倍数
nfft = fscanf(fid,'%d',1); %FFT长度
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[1,inf]); %按行读入原始信号数据
status=fclose(fid);
%计算输入数据长度
nt=length(x);
%最大细化截止频率
fa=fi+0.5*fs/np;
%取大于nt且与nt最接近的2的整数次方为FFT长度
nf=2^nextpow2(nt);
%确定细化带宽的数据长度
na=round(0.5*nf/np+1);
%频移
%建立一个按1递增的向量
n=0:nt-1;
%确定单位旋转因子向量
b=n*pi*(fi+fa)/fs;
%乘以单位旋转因子进行频移
y=x.*exp(-i*b);
%频移后的低通滤波（频域滤波）
%FFT变换
b=fft(y,nf);
%正频率带通内的元素赋值
a(1:na)=b(1:na);
%负频率带通内的元素赋值
a(nf-na+1:nf)=b(nf-na+1:nf);
%IFFT变换
b=ifft(a,nf);
%重采样
c=b(1:np:nt);
%进行细化FFT变换
a=fft(c,nfft)*2/nfft;
%变换结果重新排序
y2=zeros(1,nfft/2);
%排列负频率段的数据
y2(1:nfft/4)=a(nfft-nfft/4+1:nfft);
%排列正频率段的数据
y2(nfft/4+1:nfft/2)=a(1:nfft/4);
n=0:(nfft/2-1);
%定义细化FFT频率向量
f2=fi+n*2*(fa-fi)/nfft;
%对输入数据做FFT用来变换比较
%FFT变换
y1=fft(x,nfft)*2/nfft;
f1=n*fs/nfft;
%定义与细化FFT频率向量相同的频率范围
ni=round(fi*nfft/fs+1);
na=round(fa*nfft/fs+1);
%绘制输入时程曲线图形
subplot(2,1,1);
t=0:1/fs:(nt-1)/fs;
plot(t,x);
xlabel('时间 (s)');
ylabel('幅值');
grid on;
%在相同的频率范围下绘制幅频曲线图
nn=ni:na;
plot(f1(nn),abs(y1(nn)),':',f2,abs(y2));
xlabel('频率 (Hz)');
ylabel('幅值');
legend('普通','细化');
grid on;
%打开文件输出细化后的数据
fid=fopen(fno,'w');
for k=1:nfft/2
%每行输出3个实型数据，f为频率，y2的实和虚部
fprintf(fid,'%f %f %fn',f2(k),real(y2(k)),imag(y2(k)));
end
status=fclose(fid);
程序6-4
%三分之一倍频程处理
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%
fni=input('三分之一倍频程处理-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[1,inf]); %按行输入数据
status=fclose(fid);
%定义三分之一倍频程的中心频率
f = [1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.300 8.00];
fc =[f,10*f,100*f,1000*f,10000*f];
%中心频率与下限频率的比值
oc6=2^(1/6);
%取中心频率总的长度
nc=length(fc);
%输入数据的长度
n=length(x);
%大于并最接近n的2的幂次方长度
nfft=2^nextpow2(n);
%FFT变换
a=fft(x,nfft);
for j=1:nc
%下限频率
fl=fc(j)/oc6;
%上限频率
fu=fc(j)*oc6;
%下限频率对应的序号
nl=round(fl*nfft/sf+1);
%上限频率对应的序号
nu=round(fu*nfft/sf+1);
%如果上限频率大于折叠频率则循环中断
if fu > sf/2
m=j-1; break
end
%以每个中心频率段为通带进行带通频域滤波
b=zeros(1,nfft);
b(nl:nu)=a(nl:nu);
b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
c=ifft(b,nfft);
%计算对应每个中心频率段的有效值
yc(j)=sqrt(var(real(b(1:n))));
end
%绘制输入时程曲线图形
subplot(2,1,1);
t=0:1/sf:(n-1)/sf;
plot(t,x);
xlabel('时间 (s)');
ylabel('加速度 (g)');
grid on;
%绘制三分之一倍频程有效值图形
subplot(2,1,2);
plot(fc(1:m),yc(1:m));
xlabel('频率 (Hz)');
ylabel('有效值');
grid on;
%打开文件输出三分之一倍频程数据
fid=fopen(fno,'w');
for k=1:m
fprintf(fid,'%f %fn',fc(k),yc(k));
end
status=fclose(fid);
程序6-5
%倒频谱函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%
fni=input('倒频谱函数-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
%分析内容(1=实倒频谱,2=复倒频谱)
fun = fscanf(fid,'%d',1);
fno = fscanf(fid,'%s',1); %输出数据文件名
x = fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
n=length(x);
%建立离散时间向量
t=0:1/sf:(n-1)/sf;
%计算倒频谱函数
if fun==1
y=rceps(x); %计算实倒频谱,
else
y=cceps(x); %计算复倒频谱
end
%绘制输入时程曲线图形
subplot(2,1,1);
plot(t,x);
xlabel('时间 (s)');
ylabel('幅值');
grid on;
%绘制倒频谱时程曲线图形
subplot(2,1,2);
plot(t,y);
xlabel('时间 (s)');
ylabel('幅值');
grid on;
%打开文件输出倒频谱的数据
fid=fopen(fno,'w');
for k=1:n
fprintf(fid,'%d %fn',t(k),y(k));
end
status=fclose(fid);
程序6-6
%地震信号加速度反应谱
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%
fni=input('地震信号加速度反应谱-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
dp = fscanf(fid,'%f',1); %阻尼比
hf = fscanf(fid,'%f',1); %计算最高截止频率
df = fscanf(fid,'%f',1); %频率间隔
fno = fscanf(fid,'%s',1); %输出数据文件名
xg = fscanf(fid,'%f',[1,inf]); %输入地震波数据存成行向量
status=fclose(fid);
%计算输入数据长度
n=length(xg);
%计算反应谱长度
m=round(hf/df)+1;
%定义频率向量
f=0:df:(m-1)*df;
%定义时间向量
t=0:1/sf:(n-1)/sf;
%计算加速度反应谱
for j=0:m-1
%计算单位脉冲反应
w=2*pi*df*j;
wd=w*sqrt(1-dp*dp);
e=exp(-t*w*dp);
a=t.*wd;
s=sin(a).*((1-2*dp*dp)/(1-dp*dp));
c=cos(a).*(2*dp/sqrt(1-dp*dp));
h=wd*e.*(s+c)/sf;
%用卷积计算加速度反应谱值
r(j+1)=max(abs(conv(xg,h)));
end
%绘制输入地震波加速度曲线
subplot(2,1,1);
plot(t,xg);
xlabel('时间 (s)');
ylabel('加速度 (g)');
grid on;
%绘制地震波加速度反应谱曲线
subplot(2,1,2);
plot(f,r);
xlabel('频率 (Hz)');
ylabel('加速度 (g)');
grid on;
%打开文件输出地震波加速度反应谱数据
fid=fopen(fno,'w');
for k=1:m
%每行输出两个实型数据，f为频率，r为反应谱值
fprintf(fid,'%f %fn',f(k),r(k));
end
status=fclose(fid);
程序7-1
%生成扫频信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all
%%%%%%%%%%%%%%%%%%%%%%
fni=input('生成扫频信号-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
f0 = fscanf(fid,'%f',1); %最小截止频率
f1 = fscanf(fid,'%f',1); %最大截止频率
t0 = fscanf(fid,'%f',1); %开始时刻(s)
t1 = fscanf(fid,'%f',1); %终止时刻(s)
am = fscanf(fid,'%f',1); %扫频信号幅值
md = fscanf(fid,'%d',1); %扫频方式(1=线性,2=对数,3=指数)
fno = fscanf(fid,'%s',1); %输出数据文件名
status=fclose(fid);
%四舍五入取整求扫频信号开始时刻对应数组元素的下标
n0=round(sf*t0)+1;
%四舍五入取整求扫频信号终止时刻对应数组元素的下标
n1=round(sf*t1)+1;
%建立离散时间向量t
t=0:1/sf:(n1-1)/sf;
%按指定的扫频方式计算随时间变化频率向量
switch md
%计算扫频方式为线性时随时间变化频率向量
case 1
ft=(t-t0)./(t1-t0);
%计算扫频方式为对数时随时间变化频率向量
case 2
ft=log(t-t0)./log(t1-t0+1);
%计算扫频方式为指数时随时间变化频率向量
case 3
ft=exp((t-t0)./(t1-t0)-1);
otherwise
ft=(t-t0)./(t1-t0);
end
%计算随时间变化的园频率向量
w=2*pi*(ft.*(f1-f0)+f0);
%建立一个长度与时间长度相同元素全为0的向量
y=zeros(1,n1);
%计算扫频信号
y(n0:n1)=am*sin(w(n0:n1).*t(n0:n1));
%绘制算扫频信号随时间变化的曲线图
plot(t,y);
xlabel('时间 (s)');
ylabel('幅值');
grid on;
%打开文件输出扫频信号数据
fid=fopen(fno,'w');
for k=1:n1
%每一行输出两个实型数据，t为时间，y为扫频信号值
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序7-2
%生成拍波信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('生成拍波信号-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
fp = fscanf(fid,'%f',1); %拍波频率
np = fscanf(fid,'%d',1); %拍波个数
mp = fscanf(fid,'%d',1); %每个拍波内的周波数
ti = fscanf(fid,'%f',1); %拍波之间的时间间隔
am = fscanf(fid,'%f',1); %拍波幅值
fno = fscanf(fid,'%s',1) %输出数据文件名
status=fclose(fid);
%建立完整拍波的离散时间向量
ts=0:1/sf:(round(sf*mp/fp)-1)/sf;
%建立拍波之间的离散时间间隔向量
t0=zeros(1,round(ti*sf)+1);
%计算拍波圆频率向量
w=2*pi*fp;
%生成单个拍波的离散信号
y1=am*sin(w*ts).*sin(0.5*w*ts/mp);
%建立时间间隔和单个拍波信号的组合向量
y=[t0,y1];
%通过循环建立完整的拍波信号向量
for k=1:np-1
y=[y,t0,y1];
end
%建立与整个拍波信号长度相同的离散时间向量
t=0:1/sf:(length(y)-1)/sf;
%绘制拍波信号随时间变化的曲线图
plot(t,y);
xlabel('时间 (s)');
ylabel('幅值');
grid on;
%打开文件输出生成的拍波信号数据
fid=fopen(fno,'w');
for k=1:length(y)
%每一行输出两个实型数据，t为时间，y为拍波信号值
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序7-3
%生成白噪声信号
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('生成白噪声信号-输入数据文件名:','s');
fid=fopen(fni,'r');
sf = fscanf(fid,'%f',1); %采样频率
fi = fscanf(fid,'%f',1); %最小截止频率
fa = fscanf(fid,'%f',1); %最大截止频率
tl = fscanf(fid,'%f',1); %时间长度(秒)
am = fscanf(fid,'%f',1); %最大幅值
fno = fscanf(fid,'%s',1); %输出数据文件名
status=fclose(fid);
%计算白噪声信号数据长度n
n=round(sf*tl)+1;
%建立信号离散时间向量t
t=0:1/sf:(n-1)/sf;
%大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%四舍五入取整求最小截止频率对应数组元素的下标
ni=round(fi*nfft/sf+1);
%四舍五入取整求最大截止频率对应数组元素的下标
na=round(fa*nfft/sf+1);
%建立一个长度为nfft/2元素全为0的向量
r=zeros(1,nfft/2);
%将r的正频率带通内的元素赋值1，生成幅值谱
r(ni:na)=1;
%生成0~2pi的随机数为随机相位
a=2*pi*rand(1,nfft/2);
%将幅值谱和相位谱转换成实部和虚部
b=r.*exp(i*a)*nfft/2;
%将正负圆频率向量组合成一个向量
a=[b,b(nfft/2:-1:1)];
%IFFT变换，并取实部为生成白噪声信号
x=real(ifft(a,nfft));
%取指定长度的数据并使数据的最大值为指定的幅值
y=am*x(1:n)/max(abs(x(1:n)));
%定义自功率谱的FFT长度
nft=1024;
%计算生成白噪声信号自功率谱
y1=psd(y,nft,sf,boxcar(nft),nft/2);
%定义自功率谱的离散频率向量
f=0:sf/nft:(nft/2-1)*sf/nft;
%绘制白噪声信号随时间变化的曲线图
subplot(2,1,1);
plot(t,y);
xlabel('时间 (s)');
ylabel('幅值');
grid on;
%绘制自功率谱曲线图
subplot(2,1,2);
plot(f(1:nft/4),y1(1:nft/4));
xlabel('频率 (Hz)');
ylabel('幅值');
grid on;
%打开文件输出生成的白噪声信号数据
fid=fopen(fno,'w');
for k=1:n
%每一行输出两个实型数据，t为时间，y为白噪声信号值
fprintf(fid,'%f %fn',t(k),y(k));
end
status=fclose(fid);
程序7-4
%生成加速度人工地震波
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
%%%%%%%%%%%%%%%%%%%%%%
fid=fopen(fni,'r');
fs = fscanf(fid,'%f',1); %采样频率
tu = fscanf(fid,'%f',1); %上升时间长度
%上升时间包络线线形(1-直线、2-抛物线、3-指数曲线)
iu = fscanf(fid,'%f',1);
%上升时间包络线线形参数(只有指数曲线需要具体参数，其余均为1)
cu = fscanf(fid,'%f',1);
ta = fscanf(fid,'%f',1); %持续时间长度
td = fscanf(fid,'%f',1); %下降时间长度
%下降时间包络线线形（1-直线、2-抛物线、3-指数曲线)
id = fscanf(fid,'%f',1);
%下降时间包络线线形参数(只有抛物线,指数曲线需要具体参数，其余为1)
cd = fscanf(fid,'%f',1);
dp = fscanf(fid,'%f',1); %阻尼比值
p = fscanf(fid,'%f',1); %概率系数(一般可取p=0.85)
nn = fscanf(fid,'%d',1); %迭代次数
fno = fscanf(fid,'%s',1);%输出数据文件名
x = fscanf(fid,'%f',[2,inf]);%反应谱频率和幅值数据
tatus=fclose(fid);
%计算地震波总的时间长度
tl=tu+ta+td;
%计算生成地震波的数据长度
nt=round(fs*tl+1);
%大于并最接近nt的2的幂次方为FFT长度
nfft=2^nextpow2(nt);
%计算频率间隔（Hz）
df=fs/nfft;
%定义反应谱的离散频率向量
f=0:df:(nfft/2-1)*df;
%计算时间间隔（s）
dt=1/fs;
%定义的离散时间向量
t=0:dt:(nt-1)*dt;
%生成0~2pi的随机数为随机相位
g=rand(1,nfft/2)*2*pi;
%建立时间包络线
%建立与地震波长度相同元素全为1的向量
en=ones(1,nt);
%上升时间阶段
%确定上升时间段的长度
l=round(tu*fs)+1;
%产生上升时间段的包络线数组元素
switch iu
case 1 %直线
en(1:l)=linspace(0,1,l);
case 2 %抛物线
a=0:l-1;
en(1:l)=(a/(l-1)).^2;
case 3 %指数曲线
a=0:l-1;
en(1:l)=1-exp(-cu*a/(l-1));
end
%持续时间阶段
%确定0时刻到持续结束时刻的时间段的长度
m=round((tu+ta)*fs)+1;
%下降时间阶段
%产生下降时间段的包络线数组元素
switch id
case 1 %直线、
en(m:nt)=linspace(1,0,nt-m+1);
case 2 %抛物线
a=0:nt-m;
en(m:nt)=1-cd*(a/(nt-m)).^2;
case 3 %指数曲线
a=0:nt-m;
en(m:nt)=exp(-cd*a/(nt-m));
end
%按线性插值建立目标反应谱离散数据
%按目标反应谱的长度生成元素全为0的向量
a0=zeros(1,nfft/2);
%取目标反应谱数据的长度
n=length(x(1,:));
%四舍五入取整求反应谱最小频率对应数组元素的下标
nb=round(x(1,1)/df)+1;
%四舍五入取整求反应谱最大频率对应数组元素的下标
ne=round(x(1,n)/df)+1;
for k=1:n-1
%四舍五入取整求反应谱前一个频率数据对应数组元素的下标
l=round(x(1,k)/df)+1;
%四舍五入取整求反应谱后一个频率数据对应数组元素的下标
m=round(x(1,k+1)/df)+1;
%线性插值产生前后两个频率数据间反应谱数组元素
a0(l:m)=linspace(x(2,k),x(2,k+1),m-l+1);
end
%根据目标反应谱计算对应的近似功率谱
a1=a0;
s=zeros(1,nfft/2);
k=nb:ne;
s(k)=2*dp/pi.*(a1(k).^2)./f(k)./(-2*log(-(log(p)*pi/tl)./f(k)));
%将功率谱转换成傅立叶幅值谱
b1=sqrt(4*df*s)*nfft/2;
%定义元素全为0的反应谱传递函数矩阵,
hf=zeros(ne,nfft);
%计算加速度反应谱传递函数矩阵
for j=0:ne-1
w=2*pi*df*j;
wd=w*sqrt(1-dp*dp);
e=exp(-t.*w*dp);
a=t.*wd;
s=sin(a).*((1-2*dp*dp)/(1-dp*dp));
c=cos(a).*(2*dp/sqrt(1-dp*dp)); '
%计算加速度反应谱的脉冲响应函数向量
h=wd*e.*(s+c)/fs;
%通过FFT变换求加速度反应谱传递函数向量
hf(j+1,:)=fft(h,nfft);
end
mm=nn;
%进行生成人工地震波迭代计算
%100为最大迭代次数
for k=1:100
%将幅值谱和相位谱转换成实部和虚部
c=b1.*exp(i*g);
%将正负圆频率傅立叶谱向量组合成一个向量
d=[c,c(nfft/2:-1:1)];
%IFFT变换，并取变换结果实部为生成的地震波
e=ifft(d,nfft);
%给生成的地震波加上强度包络线
y=en.*real(e(1:nt));
%计算反应谱
%对生成的地震波进行FFT变换
yf=fft(y,nfft);
for j=1:ne
%用地震波FFT变换结果和反应谱传递函数的乘积的逆变换做卷积运算
d=ifft(yf.*hf(j,:),nfft);
%求各频率对应地震的最大响应
a1(j)=max(real(d(1:nt)));
end
%如果达到指定迭代次数显示图形
if k==mm
subplot(2,1,1);
%同时显示生成的地震波和强度包络线
plot(t,y,t,en,t,-en);
xlabel('时间 (s)');
ylabel('加速度 (g)');
grid on;
subplot(2,1,2);
%同时显示期望反应谱与反应计算谱
l=1:ne;
plot(f(l),a0(l),':',f(l),a1(l));
xlabel('频率 (Hz)');
ylabel('加速度 (g)');
legend('目标谱','计算谱');
grid on;
ig=input('继续迭代次数[取值1-9，否则退出]：');
if ig>0 & ig<10 %如果输入数字是1-9
mm=mm+ig; %增加迭代次数1-9
else %否则
break; %中断循环
end
end
c=b1;
%期望谱与计算谱的比值来修改傅立叶幅值谱
j=nb:ne;
b1(j)=c(j).*a0(j)./a1(j);
end
%打开文件输出人工地震波数据
fid=fopen(fno,'w');
for k=1:nt
%每一行输出两个实型数据，t为时间，y为人工地震波信号值
fprintf(fid,'%f %fn',t(k),y1(k));
end
status=fclose(fid);
程序8-1
%导纳圆法模态参数识别
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%
%声明为全局变量
global mn;
fni=input('导纳圆法模态参数识别-输入文件名:','s');
fid=fopen(fni,'r');
mn = fscanf(fid,'%d',1); %模态阶数
df = fscanf(fid,'%f',1); %频率间隔
f0 = fscanf(fid,'%f',[1,mn]); %模态频率初值数组
fno = fscanf(fid,'%s',1); %输出数据文件名
H0 = fscanf(fid,'%f',[2,inf]); %实测频响函数实部虚部数据
status=fclose(fid);
%定义的离散频率向量
f=0:df:(length(H0(1,:))-1)*df;
%定义的离散圆频率向量
w=2*pi*f;
%建立模态参数向量
fds=zeros(1,3*mn);
%进行4次迭代分离频响函数各模态
for n=1:4
for j=1:mn
%消除其它模态的叠加影响
l=3*(j-1);
fds(l+1:l+3)=zeros(1,3);
H1=fun81(fds,w);
H=H0-H1;
%确定拟合导纳圆的数据
k=round(f0(j)/df)+1;
[m,kc]=max(abs(H(2,k-4:k+4)));
kc=kc+k-5;
%设频响函数虚部模态峰值的0.15作为取数据的条件
v=0.15*abs(H(2,kc));
%依照条件进行取数据
for k=2:30
l=kc-k;
if abs(H(2,l))>v & abs(H(2,l))<=abs(H(2,l+1))
k1=l;
else
break;
end
end
for k=2:30
l=kc+k;
if abs(H(2,l))>v & abs(H(2,l))<=abs(H(2,l-1))
k2=l;
else
break;
end
end
%将确定的频响函数的实部和虚部数据分别存入x和y
x=H(1,k1:k2);
y=H(2,k1:k2);
%用最小二乘法进行导纳圆的拟合
A=[sum(x.^2),sum(x.*y),sum(x);
sum(x.*y),sum(y.^2),sum(y);
sum(x),sum(y),length(x)];
B=[-sum(x.^3+x.*y.^2);
-sum(x.^2.*y+y.^3);
-sum(x.^2+y.^2)];
%解线性方程
c=inv(A)*B;
%计算圆心X坐标
x0=-c(1)/2;
%计算圆心Y坐标
y0=-c(2)/2;
%计算圆半径
r=sqrt(x0^2+y0^2-c(3));
%用抛物线插值法计算模态频率值
%y1、y2、y3分别频响函数虚部模态峰值及左右的3个点
%x1、x2、x3是对应于y1、y2、y3的频率值
x1=f(kc-1);
x2=f(kc);
x3=f(kc+1);
y1=H(2,kc-1);
y2=H(2,kc);
y3=H(2,kc+1);
fc=((y2-y1)*(x3^2-x2^2)-(y3-y2)*(x2^2-x1^2))/(2*(2*y2-y1-y3)*(x2-x1));
yc=(y1*(x2-fc)^2-y2*(x1-fc)^2)/((x2-fc)^2-(x1-fc)^2);
if abs(y1)>abs(y3)
xc=H(1,kc-1)+(H(1,kc)-H(1,kc-1))*(fc-x1)/(x2-x1);
else
xc=H(1,kc)+(H(1,kc+1)-H(1,kc))*(fc-x2)/(x3-x2);
end
F(j)=fc;
%计算对应模态频率的阻尼比值
[m,k]=max(abs(y));
%根据三角形的余弦定理计算峰值点与其左右点的夹角
a=(xc-x(k-1))^2+(yc-y(k-1))^2;
b=(x0-xc)^2+(y0-yc)^2;
c=(x0-x(k-1))^2+(y0-y(k-1))^2;
a1=acos((b+c-a)/(2*sqrt(b)*sqrt(c)))/2;
a=(x(k+1)-xc)^2+(y(k+1)-yc)^2;
b=(x0-x(k+1))^2+(y0-y(k+1))^2;
c=(x0-xc)^2+(y0-yc)^2;
a2=acos((b+c-a)/(2*sqrt(b)*sqrt(c)))/2;
D(j)=(f(l+1)-f(l-1))/(F(j)*(tan(a1)+tan(a2)));
%计算对应模态频率的振型系数
S(j)=yc*2*D(j);
l=3*(j-1);
fds(l+1:l+3)=[2*pi*F(j) D(j) S(j)];
end
end
%计算拟合的频响函数
H1=fun81(fds,w);
%绘制频响函数实部拟合曲线图
figure(1)
subplot(2,1,1);
plot(f,H0(1,:),':',f,H1(1,:));
xlabel('频率 (Hz)');
ylabel('实部');
legend('实测','拟合');
grid on;
%绘制频响函数虚部拟合曲线图
subplot(2,1,2);
plot(f,H0(2,:),':',f,H1(2,:));
xlabel('频率 (Hz)');
ylabel('虚部');
legend('实测','拟合');
grid on;
%绘制频响函数导纳圆拟合曲线图
figure(2)
plot(H0(1,:),H0(2,:),':',H1(1,:),H1(2,:));
axis('equal');
xlabel('实部');
ylabel('虚部');
legend('实测','拟合');
grid on;
%打开文件输出模态参数数据
fid=fopen(fno,'w');
fprintf(fid,' 频率(Hz) 阻尼比(%%) 振型系数n');
for k=1:mn
l=(k-1)*3;
fprintf(fid,'%10.4f %10.4f %10.6fn',fds(l+1),fds(l+2)*100.0,fds(l+3));
end
status=fclose(fid);
函数8-1(用于程序8-1)
function m=fun81(x,w)
%通过实模态参数计算拟合频响函数
%输入参数
%x-实模态参数向量
%w-频率变量向量
%输出参数
%m-拟合频响函数向量
global mn
m=zeros(2,length(w));
for k=0:mn-1
l=k*3;
if x(l+1)==0
continue;
end
x1=x(l+1); %模态频率
x2=x(l+2); %阻尼比
x3=x(l+3); %振型系数
%计算拟合频响函数的实部
m(1,:)=m(1,:)-x3*(x1^2-w.^2).*w.^2./((x1^2-w.^2).^2+(2*x2*x1*w).^2);
%计算拟合频响函数的虚部
m(2,:)=m(2,:)+2*x3*x2*x1*w.^3./((x1^2-w.^2).^2+(2*x2*x1*w).^2);
end
程序8-2
%最小二乘法模态参数识别（复模态-频响函数实虚部）
%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%%
%声明全局变量
global mn
fni=input('最小二乘法模态参数识别-输入数据文件名:','s');
fid=fopen(fni,'r');
mn = fscanf(fid,'%d',1); %模态阶数
df = fscanf(fid,'%f',1); %频率间隔
f0 = fscanf(fid,'%f',mn); %输入模态频率初值数组
d0 = fscanf(fid,'%f',mn); %输入模态阻尼比初值数组
fno = fscanf(fid,'%s',1); %输出数据文件名
b = fscanf(fid,'%f',[2,inf]); %实测频响函数实部虚部数据
status=fclose(fid);
%建立离散频率向量
f=0:df:(length(b(1,:))-1)*df;
%建立离散圆频率向量
w=2*pi*f;
%建立实测频响函数复数向量
H =b(1,:)+b(2,:)*i;
%计算模态圆频率初值向量
w0=2*pi*f0;
%建立模态初参数向量
for j=1:mn
l=4*(j-1);
x0(l+1:l+4)=[-w0(j)*d0(j),w0(j)*sqrt(1-d0(j)^2),1,1];
end
%用最小二乘非线性数据拟合法估计复模态参数
x=lsqcurvefit('fun82',x0,w,H);
%计算模态频率、阻尼比和振型系数
for j=1:mn
l=4*(j-1);
c=x(l+1)+i*x(l+2);
d=x(l+3)+i*x(l+4);
F(j)=abs(c)/(2*pi); %模态频率
D(j)=-real(c)/abs(c);%阻尼比
S(j)=d; %复振型系数
end
%计算拟合的频响函数向量
H1=fun82(x,w);
%绘制频响函数实部拟合曲线图
subplot(2,1,1);
plot(f,real(H),':',f,real(H1));
xlabel('频率 (Hz)');
ylabel('实部');
legend('实测','拟合');
grid on;
%绘制频响函数虚部拟合曲线图
subplot(2,1,2);
plot(f,imag(H),':',f,imag(H1),'r');
xlabel('频率 (Hz)');
ylabel('虚部');
legend('实测','拟合');
grid on;
%打开文件输出模态参数数据
fid=fopen(fno,'w');
fprintf(fid,' 频率(Hz) 阻尼比(%%) 振型系数n');
for k=1:mn
fprintf(fid,'%10.4f %10.4f %10.6fn',F(k),D(k)*100.0,imag(S(k)));
end
status=fclose(fid);
函数8-2(用于程序8-2)
function m=fun82(x,w)
%通过模态参数计算拟合频响函数值
%输入参数
%x-复模态参数向量
%w-频率变量向量
%输出参数
%m-拟合频响函数向量
global mn
m=zeros(1,length(w));
for k=0:mn-1
l=4*k;
x1=x(l+1); %特征值实部
x2=x(l+2); %特征值虚部
x3=x(l+3); %特征向量实部
x4=x(l+4); %特征向量虚部
m=m-w.^2.*((x3+i*x4)./(w*i-(x1+i*x2))+(x3-i*x4)./(w*i-(x1-i*x2)));
end
程序8-2a
%最小二乘法模态参数识别（实模态-频响函数虚部）
%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%
%声明为全局变量
global mn
fni=input('最小二乘法模态参数识别-输入数据批文件名:','s');
fid=fopen(fni,'r');
mn = fscanf(fid,'%d',1); %模态阶数
df = fscanf(fid,'%f',1); %频率间隔
f0 = fscanf(fid,'%f',mn); %输入模态频率初值数组
d0 = fscanf(fid,'%f',mn); %输入模态阻尼比初值数组
s0 = fscanf(fid,'%f',mn); %输入模态振型系数初值数组
fno = fscanf(fid,'%s',1); %输出数据文件名
b = fscanf(fid,'%f',[2,inf]); %输入实测频响函数实虚部数据数组
status=fclose(fid);
%取频响函数长度
n=length(b(1,:));
%建立离散频率向量
f=0:df:(n-1)*df;
%建立离散圆频率向量
w=2*pi*f;
%建立实测频响函数虚部向量
H =b(2,:);
%建立模态圆频率初值向量
w0=2*pi*f0;
%建立模态初参数向量
for j=1:mn
l=3*(j-1);
x0(l+1:l+3)=[w0(j) d0(j) s0(j)];
end
%用最小二乘非线性拟合法估计实模态参数
x=lsqcurvefit('fun82a',x0,w,H);
%计算模态频率、阻尼比和振型系数
for j=1:mn
l=3*(j-1);
F(j)=x(l+1)/(2*pi);%模态频率
D(j)=x(l+2); %阻尼比
S(j)=x(l+3); %振型系数
end
%建立实测频响函数复数向量
H=b(1,:)+i*b(2,:);
%计算拟合频响函数的实部
c=zeros(1,n);
for k=0:mn-1
l=k*3;
x1=x(l+1);
x2=x(l+2);
x3=x(l+3);
c=c+x3*(x1^2-w.^2).*w.^2./((x1^2-w.^2).^2+(2*x1*x2*w).^2);
end
%计算拟合频响函数的虚部
d=fun82a(x,w);
%建立拟合频响函数复数向量
H1=c+i*d;
%绘制频响函数实部拟合曲线图
subplot(2,1,1);
plot(f,real(H),':',f,real(H1));
xlabel('频率 (Hz)');
ylabel('实部');
legend('实测','拟合');
grid on;
%绘制频响函数虚部拟合曲线图
subplot(2,1,2);
plot(f,imag(H),':',f,imag(H1),'r');
xlabel('频率 (Hz)');
ylabel('虚部');
legend('实测','拟合');
grid on;
%打开文件输出模态参数数据
fid=fopen(fno,'w');
fprintf(fid,' 频率(Hz) 阻尼比(%%) 振型系数n');
for k=1:mn
fprintf(fid,'%10.4f %10.4f %10.6fn',F(k),D(k)*100.0,S(k));
end
status=fclose(fid);
函数8-2a(用于程序8-2a)
function m=funlsq(x,w)
%通过实模态参数计算拟合频响函数虚部
%输入参数
%x-实模态参数向量
%w-频率变量向量
%输出参数
%m-拟合频响函数虚部向量
global mn
m=zeros(1,length(w));
for k=0:mn-1
l=k*3;
x1=x(l+1); %模态频率
x2=x(l+2); %阻尼比
x3=x(l+3); %振型系数
m=m-2*x3*x2*x1*w.^3./((x1^2-w.^2).^2+(2*x2*x1*w).^2);
end
程序8-3
%加权最小二乘法模态参数识别（复模态-频响函数实虚部）
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%
%声明为全局变量
global mn w H W ig
fni=input('加权最小二乘法模态参数识别-输入数据文件名:','s');
fid=fopen(fni,'r');
mn = fscanf(fid,'%d',1); %模态阶数
df = fscanf(fid,'%f',1); %频率间隔
f0 = fscanf(fid,'%f',mn); %输入模态频率初值数组
d0 = fscanf(fid,'%f',mn); %输入模态阻尼比初值数组
fno = fscanf(fid,'%s',1); %输出数据文件名
b = fscanf(fid,'%f',[2,inf]);%输入实测频响函数实部虚部数据数组
status=fclose(fid);
n=length(b(1,:));
%建立离散频率向量
f=0:df:(n-1)*df;
%建立离散圆频率向量
w=2*pi*f;
%建立实测频响函数复数向量
H =b(1,:)+b(2,:)*i;
%计算模态圆频率初值向量
w0=2*pi*f0;
%建立模态初参数向量
for j=1:mn
l=4*(j-1);
x0(l+1:l+4)=[-w0(j)*d0(j),w0(j)*sqrt(1-d0(j)^2),1,1];
end
%建立权向量
W=(abs(H)./max(abs(H))).^0.2;
%%用加权非线性最小二乘求解法估计复模态参数
ig=1;
x=lsqnonlin('fun83',x0);
%计算模态频率、阻尼比和振型系数
for j=1:mn
l=4*(j-1);
c=x(l+1)+i*x(l+2);
d=x(l+3)+i*x(l+4);
F(j)=abs(c)/(2*pi); %模态频率
D(j)=-real(c)/abs(c);%阻尼比
S(j)=d; %复振型系数
end
%计算拟合频响函数
ig=0;
H1=fun83(x);
%绘制频响函数实虚部拟合曲线图
subplot(2,1,1);
plot(f,real(H),':',f,real(H1));
xlabel('频率 (Hz)');
ylabel('实部');
legend('实测','拟合');
grid on;
subplot(2,1,2);
plot(f,imag(H),':',f,imag(H1));
xlabel('频率 (Hz)');
ylabel('虚部');
legend('实测','拟合');
grid on;
%打开文件输出模态参数数据
fid=fopen(fno,'w');
fprintf(fid,' 频率(Hz) 阻尼比(%%) 振型系数n');
for k=1:mn
fprintf(fid,'%10.4f %10.4f %10.6fn',F(k),D(k)*100.0,imag(S(k)));
end
status=fclose(fid);
函数8-3(用于程序8-3)
function m=fun83(x)
%通过复模态参数计算拟合频响函数
%x-复模态参数向量
%w-频率变量向量
%H-实测频响函数向量
%W-权函数向量
%ig=0计算拟合频响函数
%ig=1计算拟合频响函数与实测频响函数的差值
global mn w H W ig
m=zeros(1,length(w));
for k=0:mn-1
l=4*k;
x1=x(l+1); %特征值实部
x2=x(l+2); %特征值虚部
x3=x(l+3); %特征向量实部
x4=x(l+4); %特征向量虚部
m=m-w.^2.*((x3+i*x4)./(w*i-(x1+i*x2))+(x3-i*x4)./(w*i-(x1-i*x2)));
end
if ig>0
m=(m-H).*W;
end
程序8-3a
%加权最小二乘法模态参数识别（实模态-频响函数虚部）
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%
%声明为全局变量
global mn w H W ig
fni=input('加权最小二乘法模态参数识别-输入数据文件名:','s');
fid=fopen(fni,'r');
mn = fscanf(fid,'%d',1); %模态阶数
df = fscanf(fid,'%f',1); %频率间隔
f0 = fscanf(fid,'%f',mn); %输入模态频率初值数组
d0 = fscanf(fid,'%f',mn); %输入模态阻尼比初值数组
s0 = fscanf(fid,'%f',mn); %输入模态振型系数初值数组
fno = fscanf(fid,'%s',1); %输出数据文件名
b = fscanf(fid,'%f',[2,inf]); %输入实测频响函数实虚部数据
status=fclose(fid);
%取频响函数长度
n=length(b(1,:));
%建立离散频率向量
f=0:df:(n-1)*df;
%建立离散圆频率向量
w=2*pi*f;
%建立实测频响函数虚部向量
H =b(2,:);
%建立模态圆频率初值向量
w0=2*pi*f0;
%建立模态初参数向量
for j=1:mn
l=3*(j-1);
x0(l+1:l+3)=[w0(j) d0(j) s0(j)];
end
%建立权向量
W=(abs(H)./max(abs(H))).^0.3;
%用加权非线性最小二乘求解法估计实模态参数
ig=1;
x=lsqnonlin('fun83a',x0);
%计算模态频率、阻尼比和振型系数
for j=1:mn
l=3*(j-1);
F(j)=x(l+1)/(2*pi); %模态频率
D(j)=x(l+2); %阻尼比
S(j)=x(l+3); %振型系数
end
%建立实测频响函数复数向量
H=b(1,:)+i*b(2,:);
%计算拟合频响函数的实部
c=zeros(1,n);
for k=0:mn-1
l=k*3;
x1=x(l+1);
x2=x(l+2);
x3=x(l+3);
c=c+x3*(x1^2-w.^2).*w.^2./((x1^2-w.^2).^2+(2*x2*x1*w).^2);
end
%计算拟合频响函数的虚部
ig=0;
d=fun83a(x);
%计算拟合频响函数复数向量
H1=c+i*d;
%绘制频响函数实部拟合曲线图
subplot(2,1,1);
plot(f,real(H),':',f,real(H1));
xlabel('频率 (Hz)');
ylabel('实部');
legend('实测','拟合');
grid on;
%绘制频响函数虚部拟合曲线图
subplot(2,1,2);
plot(f,imag(H),':',f,imag(H1));
xlabel('频率 (Hz)');
ylabel('虚部');
legend('实测','拟合');
grid on;
%打开文件输出模态参数数据
fid=fopen(fno,'w');
fprintf(fid,' 频率(Hz) 阻尼比(%%) 振型系数n');
for k=1:mn
fprintf(fid,'%10.4f %10.4f %10.6fn',F(k),D(k)*100.0,S(k));
end
status=fclose(fid);
函数8-3a(用于程序8-3a)
function m=fun83(x)
%通过实模态参数计算拟合频响函数虚部及与实测值的差值
%输入参数
%x-实模态参数向量
%w-频率变量向量
%ig=0计算拟合频响函数虚部
%ig=1计算拟合频响函数与实测频响函数虚部的差值
%输出参数
%拟合频响函数向量或差值向量
global mn w H W ig
m=zeros(1,length(w));
for k=0:mn-1
l=k*3;
x1=x(l+1); %模态频率
x2=x(l+2); %阻尼比
x3=x(l+3); %振型系数
m=m-2*x3*x2*x1*w.^3./((x1^2-w.^2).^2+(2*x2*x1*w).^2);
end
if ig>0
m=(m-H).*W;
end
程序8-4
%有理分式多项式法模态参数识别
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('有理分式多项式法模态参数识别-输入数据文件名:','s');
fid=fopen(fni,'r');
mn = fscanf(fid,'%d',1); %模态阶数
df = fscanf(fid,'%f',1); %频率间隔
fno = fscanf(fid,'%s',1); %输出数据文件名
b = fscanf(fid,'%f',[2,inf]); %实测频响函数实、虚部数据
status=fclose(fid);
%定义幂多项式的阶次
nm = mn*2;
%取频响函数的长度
n=length(b(1,:));
%建立离散频率向量
f=0:df:(n-1)*df;
%建立离散圆频率向量
w=2*pi*f;
%建立归一化离散频率向量
wi = w/max(w);
%建立实测频响函数复数向量
H=b(1,:)+i*b(2,:);
%计算拟合频响函数的分子和分母系数向量
[A,B] = invfreqs(H,wi,nm,nm,[],100);
%幂多项式方程求根(零点)
P = roots(B);
%计算模态频率向量
F1 = abs(P)*max(w)/(2*pi);
%计算阻尼比向量
D1 = -real(P)./(abs(P));
%计算振型系数向量
for k=1:nm
if k==1
p(1:nm-1)=P(2:nm);
else
p(1:k-1)=P(1:k-1);
p(k:nm-1)=P(k+1:nm);
end
y=poly(p);
S1(k)=polyval(A,P(k))/polyval(y,P(k));
end
%计算拟合的频响函数
H1 = freqs(A,B,wi);
%绘制频响函数实部拟合曲线图
figure(2)
nn=1:n;
subplot(2,1,1);
plot(f,real(H(nn)),':',f,real(H1(nn)),'r');
xlabel('频率 (Hz)');
ylabel('实部');
legend('实测','拟合');
grid on;
%绘制频响函数虚部拟合曲线图
subplot(2,1,2);
plot(f,b(2,:),':',f,imag(H1(nn)),'r');
xlabel('频率 (Hz)');
ylabel('虚部');
legend('实测','拟合');
grid on;
%将模态频率从小到大排序
[F2,I]=sort(F1);
%剔除方程解中的非模态项(非共轭根)和共轭项(重复项)
m=0;
for k=1:nm-1
if F2(k)~=F2(k+1)
continue;
end
m=m+1;
l=I(k);
F(m)=F1(l); %模态频率
D(m)=D1(l); %阻尼比
S(m)=S1(l); %振型系数
end
%打开文件输出模态参数数据
fid=fopen(fno,'w');
fprintf(fid,' 频率(Hz) 阻尼比(%%) 振型系数n');
for k=1:m
fprintf(fid,'%10.4f %10.4f %10.6fn',F(k),D(k)*100.0,imag(S(k)));
end
status=fclose(fid);
程序8-5
%正交多项式法模态参数识别
%%%%%%%%%%%%%%%%%%%%%%%%%
clear
clc
close all hidden
format long
%%%%%%%%%%%%%%%%%%%%%%%%%
fni=input('正交多项式法模态参数识别-输入数据文件名:','s');
fid=fopen(fni,'r');
mn = fscanf(fid,'%d',1); %模态阶数
df = fscanf(fid,'%f',1); %频率间隔
fno = fscanf(fid,'%s',1); %输出数据文件名
h = fscanf(fid,'%f',[2,inf]); %实测频响函数实部虚部数据
status=fclose(fid);
%定义正交多项式的阶次
nm=2*mn;
%取频响函数的长度
n=length(h(1,:));
%建立离散频率向量
f=0:df:(n-1)*df;
%建立从负到正的归一化频率向量
w=[-f(n:-1:1) f(1:n)]/max(f);
%建立对应频率向量的实测频响函数向量
H=[h(1,n:-1:1)-i*h(2,n:-1:1) h(1,1:n)+i*h(2,1:n)];
%将频率向量设置为列向量
[k,l]=size(w);
if k w=w.'; end %将频响函数向量设置为列向量 [k,l]=size(H); if k H=H.'; end %建立分子正交多项式及幂多项式的转换矩阵 [p,cp]=fun85(H,w,1,nm); %建立分母正交多项式及幂多项式的转换矩阵 [q,cq]=fun85(H,w,2,nm); %最小二乘法求解正交多项式的分子系数C和分母向量D P=p; for k=1:nm Q(:,k)=(H.*q(:,k)); end T=-real(P'*Q); g=real(P'*(H.*q(:,nm+1))); D=-inv(eye(nm)-T'*T)*T'*g; C=g-T*D; D(nm+1)=1; %通过转换矩阵计算分子幂多项式系数向量A A=cp*C; %将A设置成按降幂排序的列向量 A=A(nm+1:-1:1).'; %通过转换矩阵计算分母幂多项式系数向量A B=cq*D; %将B转换成按降幂排序的列向量 B=B(nm+1:-1:1).'; %计算有理分式多项式的极点和留数 [R,U,K]=residue(A,B); %计算模态频率向量 F1 = abs(U)*max(f); %计算阻尼比向量 D1 = -real(U)./(abs(U)); %计算振型系数向量 for k=1:nm if k==1 u(1:nm-1)=U(2:nm); else u(1:k-1)=U(1:k-1); u(k:nm-1)=U(k+1:nm); end y=poly(u); S1(k)=polyval(A,U(k))/polyval(y,U(k)); end %建立正频率的频响函数向量 H=h(1,1:n)+i*h(2,1:n); %计算生成的频响函数 for l=1:n H1(l)=sum(C'.*p(n+l,:))/sum(D'.*q(n+l,:)); end %绘制频响函数实部拟合曲线图 subplot(2,1,1); plot(f,real(H),':',f,real(H1),'r'); xlabel('频率 (Hz)'); ylabel('实部'); legend('实测','拟合'); grid on; %绘制频响函数虚部拟合曲线图 subplot(2,1,2); plot(f,imag(H),':',f,imag(H1),'r'); xlabel('频率 (Hz)'); ylabel('虚部'); legend('实测','拟合'); grid on; %将模态频率从小到大排序 [F2,I]=sort(F1); %剔除方程解中的非模态项(非共轭根)和共轭项(重复项) m=0; for k=1:nm-1 if F2(k)~=F2(k+1) continue; end m=m+1; l=I(k); F(m)=F1(l); %模态频率 D(m)=D1(l); %阻尼比 S(m)=S1(l); %振型系数 end %打开文件输出模态参数数据 fid=fopen(fno,'w'); fprintf(fid,' 频率(Hz) 阻尼比(%%) 振型系数n'); for k=1:m fprintf(fid,'%10.4f %10.4f %10.6fn',F(k),D(k)*100.0,imag(S(k))); end status=fclose(fid); 函数8-5(用于程序8-5) function [P,cf]=fun85(H,w,ig,m) %建立分子或分母正交多项式矩阵及相应幂多项式的转换矩阵 %输入参数 % H = 频响函数向量 % w = 频率向量 % ig = 1为计算分子多项式，=2为计算分母多项式 % m = 多项式的阶次 %输出参数 % P = 正交多项式矩阵 % cf = 转换正交多项式系数成幂多项式系数的转换矩阵 % %建立分子或分母正交多项式矩阵的权向量 if ig==1 %计算分子正交多项式矩阵的权向量 W=ones(size(w)); else %计算分母正交多项式矩阵的权向量 W=(abs(H)).^2;
版权声明：本文标题：王济matlab在振动信号处理中的应用代码 内容由网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://www.roclinux.cn/b/1709776758a546229.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论