跪求一个matlab 求取瞬时振幅频率和相位的程序么，输入一个信号，输出这三个，谢谢

我最近也在做这个，我用的是1kHz单频正弦波

unwrap的作用是解卷绕

没用unwrap函数时，瞬时相位是这样的

加上unwrap函数之后是这样的

fs=80是你信号的采样频率

xhd1=fs*diff(xh1)/(2*pi)是要求信号的瞬时频率，就是瞬时相位的导数

前面乘fs是为了跟实际频率对应上

不乘fs是这样的

纵轴频率对应的不是1kHz，但是乘了fs就变成下面这样，纵轴频率就是实际频率了

想换别的信号，就把xn换成你想改的信号，再根据你的采样频率改一下fs就行了

希尔伯特变换的意义本文不提，本文的目标是举例说明到底离散傅里叶变换是如何实现的，并编写程序与matlab自带的hilbert函数结果对比，验证我们的实现过程是否正确。

原始信号为 ，其离散希尔伯特变换的定义公式为：

说明：先让原始信号与 信号做卷积，然后一起合并成一个 复数 信号 。

问题：看上去很简单，但这里的卷积不是一般意义上的卷积的 *** 作！

所以：实际中得到 的方法是通过借助 离散傅里叶变换DFT 来实现的！

因此：本文就用 离散傅里叶变换 来实现 离散希尔伯特变换 。

设原始信号为 ，总长度一般 (下标n是从 1 开始到N)，总体实现步骤可分为3步：

相应的matlab程序：

结果：

手动实现正确！

希尔伯特变换的结果是给原始信号 提供了一个幅值、频率不变，但相位平移90°的信号 。

所以，希尔伯特变换是从"时域"到"时域"的变换！只改变了相位，所以又叫90°移相滤波器；

所以，原始信号 与它的希尔伯特变换 构成正交副。

现在，我们换回到最初的记法：原始信号和它对应的希尔伯特变换信号分别用 和 表示，那么对应的" 解析信号 "就可以用这两个东西组成：

对于这个解析信号，我们可以得它的3瞬属性：瞬时振幅、瞬时相位、瞬时频率。

可以看出，3瞬属性是相互关联的！

瞬时振幅、瞬时相位可以直接求，有意义；

但是瞬时频率 直接 根据解析信号这么按公式 求 ，是 没有物理意义 的！

并且 离散信号 ，它的瞬时频率求导只能按照" 差分 "来近似，即：

给出matlab的实现程序：

效果：

3瞬属性中的瞬时频率，很明显可以看出它有很多的" 负频率 "！这很明显是错误的。

所以，直接根据" 解析信号 "算瞬时频率是无意义的！

所以，真正做 3瞬属性 的分析，做原信号的" 时频谱 "分析，我们用的：

—— 希尔伯特-黄变换。

代码的解释（供参考，下列是根据经验大概的判断，如不对请谅解）

%自定义函数 spacemodel（）作用是当flag某值，执行相当于的程序命令

function [sys,x0,str,ts] = spacemodel(t,x,u,flag)

switch flag,

case 0,

[sys,x0,str,ts]=mdlInitializeSizes

case 1,

sys=mdlDerivatives(t,x,u)

case 3,

sys=mdlOutputs(t,x,u)

case {2,4,9}

sys=[]

otherwise

error(['Unhandled flag = ',num2str(flag)])

end

%自定义函数 mdlInitializeSizes（）作用是当flag=0，执行相当于的程序命令

function [sys,x0,str,ts]=mdlInitializeSizes

sizes = simsizes

sizes.NumContStates = 4

sizes.NumDiscStates = 0圆盘

sizes.NumOutputs = 4输出

sizes.NumInputs = 2输入

sizes.DirFeedthrough = 0反馈

sizes.NumSampleTimes = 1采样时间

sys = simsizes(sizes)大小

x0 = [3001] 初值

str = []空矩阵

ts = [0 0] 初值

%自定义函数 mdlDerivatives（）作用是（flag=1）根据t,x,u的值，将t,x,u变量值赋值给相对应的sys变量

function sys=mdlDerivatives(t,x,u)

%a=1

a=1000*rands(1)

d1=a*0.3*sin(3*pi*t)随机作用力

d2=a*0.1*(1-exp(-t))随机作用力

dq1=x(1)%Angle1 speed:dq1 角速度

dq2=x(2)%Angle2 speed:dq2 角速度

q1=x(3)%Angle1:q1 角度

q2=x(4)%Angle2:q2 角度

tol1=u(1)

tol2=u(2)

%%%%%%%%%%%%%%%%%%%%%%%%% From aslpd_con_s.m

m1=10m2=5 %质量

l1=1l2=0.5%长度

r1=0.5r2=0.25%半径

i1=0.83+m1*r1^2+m2*l1^2%惯性半径

i2=0.3+m2*r2^2% %惯性半径

g=9.8%重力加速度

D=[i1+i2+2*m2*r2*l1*cos(q2) i2+m2*r2*l1*cos(q2)i2+m2*r2*l1*cos(q2) i2]质量

C=[-m2*r2*l1*dq2*sin(q2),-m2*r2*l1*(dq1+dq2)*sin(q2)m2*r2*l1*dq1*sin(q2),0]阻尼

G=[(m1*r1+m2*l1)*g*cos(q1)+m2*r2*g*cos(q1+q2)m2*r2*g*cos(q1+q2)]重量

%%%%%%%%%%%%%%%%%%%%%%%%%%

D2=inv(D)D的逆矩阵

Ta=[d1d2]

A=-D2*C振幅

Z=-D2*G相位角

sys(1)=A(1,1)*x(1)+A(1,2)*x(2)+Z(1)+D2(1,1)*(-Ta(1)+tol1)+D2(1,2)*(-Ta(2)+tol2)动载荷

sys(2)=A(2,1)*x(1)+A(2,2)*x(2)+Z(2)+D2(2,1)*(-Ta(1)+tol1)+D2(2,2)*(-Ta(2)+tol2)动载荷

sys(3)=x(1) 角速度

sys(4)=x(2) 角速度

%自定义函数 mdlOutputs（）作用是（flag=3）根据t,x,u的值，将x变量值赋值给相对应的sys变量

function sys=mdlOutputs(t,x,u)

sys(1)=x(1)%Angle1 speed:dq1 角速度

sys(2)=x(2)%Angle2 speed:dq2 角速度

sys(3)=x(3)%Angle1:q1 角度

sys(4)=x(4)%Angle2:q2 角度

---