正弦扫频信号如何提取每个对应频率的幅值

三水觞

有这样一个响应信号，我要得到每个频率对应的幅值，如何求得，
下面是matlab代码
clear all;
%clf;
clc;
%%%%%%%%%梁相关参数求得频响
L=1.4;   %梁的长度
b=0.03;  %梁的宽度
h=0.012;   %梁的高度
p=7.8*10^3;   %梁所选材料为钢时的密度
E=2.1*10^11;    % 梁所选材料为钢时的弹性模量
A=b*h;     %梁的截面积
I=b*h^3/12;    %极惯性矩
kesi=0.01;
dof=10;
figure
%%%%%%%倍数扫频相关参数计算
fl=20;
fm=2000;
oct=log2(fm/fl);
T=200;
n=120000;
dt=T/n;
N=log(fm/fl)/log(2);
R=N/T;
t=0:dt:200;
fpu=fl*2.^(R*t);
fi=2*pi*fl*(-1+2.^(R*t))/(R*log(2));
% y=sin(fi);
% plot(t,y);
c=2*dof-1;
d=2*dof-11;
[K,M,C,~,Mr,~,~,~,~,~,~,~,ha19] = KMC(L,b,h,p,E,kesi,dof,c,c,fpu);
Ha19=ha19;
R1=2;
a=R1./Ha19;
f=abs(a).*sin(fi+angle(Ha19));
plot(t,f);
A=[zeros(2*dof),eye(2*dof);-inv(M)*K,-inv(M)*C];
B0=zeros(2*dof,1);
B0(c,1)=1;
B=[zeros(2*dof);inv(M)]*B0;
sys=ss(A,B,A,B);
Ftt=zeros(n+1,1);
Ftt(1:n+1,1)=Mr(c,c)*f;
x0=zeros(4*dof,1);
[y,T] = lsim(sys,Ftt,t,x0);
a1=y(:,2*dof+c);
at19=a1.';
figure
plot(t,at19)



KMc如下
%%%%%%%%%%%求解刚度矩阵，质量阵，阻尼阵，固有频率，以及位移及速度频响函数
function [K,M,C,Kr,Mr,Cr,wn,wd,fi,fh,fh1,Hwb,Ha] = KMC( L,b,h,p,E,kesi,dof,c,d,omega)
%L:梁的长度
%b:梁的宽度
%h:%梁的高度
%p:%梁所选材料为钢时的密度
%E:梁所选材料为钢时的弹性模量
%kesi:阻尼比
%dof:将梁分成的小段数目

A=b*h;     %梁的截面积
I=b*h^3/12;    %极惯性矩

%解析解
f=@(x) cos(x)*cosh(x)+1; %define the function 固有频率方程
i=0;
j=1;    %To get the lower n roots of f(x)=0.
while j<=dof
if f(i)*f(i+1)<= 0
beta(j)=fzero(f,i)  ; %cosx*coshx+1=0的解
fh(j)=beta(j)^2*sqrt(E*I/(p*A*L^4))/(2*pi);  %解析解
w=fh.*2*pi;
j=j+1;
end
i=i+1;
end
l=L/dof;       %梁分成dof段后每个小梁的长度
%仿真解
K=zeros(2*(dof+1),2*(dof+1));%总体刚度阵
M=zeros(2*(dof+1),2*(dof+1));%总体质量阵
Ke=E*I/l^3*[12 6*l -12 6*l;6*l 4*l^2 -6*l 2*l^2;
    -12 -6*l 12 -6*l;6*l 2*l^2 -6*l 4*l^2];  %每个单元梁的刚度矩阵
Me=p*A*l/420*[156 22*l 54 -13*l;22*l 4*l^2 13*l -3*l^2;
    54 13*l 156 -22*l;-13*l -3*l^2 -22*l 4*l^2];  %每个单元量的质量矩阵
for a=1:dof
    for i=1:4
        for j=1:4  
            K(i+(a-1)*2,j+(a-1)*2)=K(i+(a-1)*2,j+(a-1)*2)+Ke(i,j);%总体刚度矩阵的拼装
            M(i+(a-1)*2,j+(a-1)*2)=M(i+(a-1)*2,j+(a-1)*2)+Me(i,j);%总体质量矩阵的拼装
        end
    end   
end
K=K(3:2*(dof+1),3:2*(dof+1));%约束条件，第一点被约束,得到刚度矩阵
M=M(3:2*(dof+1),3:2*(dof+1));%约束条件，第一点被约束，得到质量矩阵
[fi,v]=eig(K,M);   %求出振型和频率
m=size(v);  % 频率矩阵的维数
Mr=fi'*M*fi; %对质量阵作模态变换   
Kr=fi'*K*fi; %对质量矩阵作模态变换   
for i=1:1:m
    wn(i)=sqrt(v(i,i));   %将频率放到一个数组中
    wd(i)=wn(i)*sqrt(1-kesi^2);%自然频率
    Cr(i)=2*kesi*wn(i)*Mr(i,i);%求Cr
end   
fh1=wn./(2*pi);  %matlab 有限元法解\
Cr=diag(Cr);%模态变换后的阻尼矩阵
%C=inv(fi')*Cr*inv(fi);
C=M*fi*inv(Mr)*Cr*inv(Mr)*fi'*M;%阻尼矩阵
%频响函数
Hw=zeros(m);
e=sqrt(-1);
Hwb=Hw(c,d);
for r=1:1:m
    Hwb=Hwb+fi(c,r)*fi(d,r)./(Kr(r,r)-(omega.*2*pi).^2.*Mr(r,r)+e*(omega.*2*pi)*Cr(r,r));
    %   可求得频响函数矩阵
end
Ha=-(omega.*2*pi).^2.*Hwb;

%%画图
%am1=subs(Hwb);
%am=abs(am1);
%plot(omega,am);
%hold on;
   %ang=angle(Ha)/pi*180;  求角度
   %amdb=10*log10(am);    频率转化为db
%xlabel('Frequence(Hz)');
%ylabel('Amp(Hz)');
%gtext('悬臂梁幅频函数元素H93(f)');%元素编号要改
%set(gcf,'color','white');%设置背景颜色为白色
%set(gca,'linewidth',1.5);%设置坐标线宽为1.5
%set(get(gca,'Children'),'linewidth',1.5);%设置线宽为1.5

end



补充内容 (2014-4-4 13:34):
简单的说就是上面图中的响应信号，如何求得它在频率内各频率的幅值

yghit08

为什么愿意贴一大串代码？
描述自己的问题就好了。想说的是从帖子题目和内容上不知问题是啥

猫头鹰先生

粤语残片

猫头鹰先生 发表于 2014-4-3 19:43

猫头鹰先生，能不能请教你http://forum.vibunion.com/forum. ... =page%3D1#pid753829这里的问题，谢谢~