符号和数值计算瞬时频率的差别

shuihai707

正常数值计算瞬时频率
fs=1000;%采样频率
N=1000;%采样点数
n=0:N-1;%采样序号
t=n/fs;%采样时间序列
x=cos(16*pi*t+12*pi*t.*t);%信号
y1=fs*abs(diff(acos(x)))/(2*pi);%求瞬时频率
plot((1:length(y1))/fs,y1);axis([0 1 0 30]);
符号计算瞬时频率
syms x t;
x=cos(12*pi*t^2+16*pi*t);
y=abs(diff(x,'t'))/(2*pi*(1-x^2)^0.5);%这个表达式和下面yy表达式是一致的，我发现通过
函数simplify(yy)求出的简化表达式是4*abs(3*t+2),所以瞬时频率就是一条直线了
% yy=abs(diff(acos(x),'t'))/(2*pi);4*abs(3*t + 2);
t0=0;
for i=1:1000
    y2(i)=subs(y,'t',t0);
    t0=t0+0.001;
end
plot((1:length(y2))/fs,y2);axis([0 1 0 30]);

高手指点迷津！！！

符号计算瞬时频率

yghit08

一个是连续的方法，一个是离散的方法

shuihai707

yghit08 发表于 2013-4-23 15:39

                               
登录/注册后可看大图

一个是连续的方法，一个是离散的方法

符号法给定的时间值是离散的啊？对于y的表达式在极值点处应该出现间断点啊？能详细说明吗？

yghit08

shuihai707 发表于 2013-4-23 16:01

                               
登录/注册后可看大图

符号法给定的时间值是离散的啊？对于y的表达式在极值点处应该出现间断点啊？能详细说明吗？

自己考虑一下。符号计算中，t不给离散数据你怎么在计算机中表达出来？利用符号计算的过程是连续的，亦即它的表达式是连续的，那么你想再做出来数值解，在计算机里就是离散的。

shuihai707

原来这样，但对于公式，图形应该出现间断点的

shuihai707

yghit08 发表于 2013-4-23 16:11

                               
登录/注册后可看大图

自己考虑一下。符号计算中，t不给离散数据你怎么在计算机中表达出来？利用符号计算的过程是连续的，亦即它 ...

原来这样，但对于公式，图形应该在极值点处出现间断点的（帮我把重复的帖子删了，操作失误，自己删不了）

yghit08

shuihai707 发表于 2013-4-23 16:39

                               
登录/注册后可看大图

原来这样，但对于公式，图形应该在极值点处出现间断点的（帮我把重复的帖子删了，操作失误，自己删不了） ...

至于间断点，你自己分析整个过程中x能否取到1.这和你数值解的情况一样，如果数值解没有出现，那么你怎么能让理论解出现间断呢？？

shuihai707

数值解很容易取到1，从数值解的图形上的突变就可以看出，突变主要就是由于数值解取到1造成的。上述x用matlab计算数值解，其实取到1的极值点只有两个，对于其他的极值点我在程序中也处理了，令极大值为1，极小值为-1（任达千论文可证，上述程序这句调试代码没写），所以对于符号图形没有出现间断点甚是疑惑。

shuihai707

yghit08 发表于 2013-4-23 17:00

                               
登录/注册后可看大图

至于间断点，你自己分析整个过程中x能否取到1.这和你数值解的情况一样，如果数值解没有出现，那么你怎么能 ...

数值解很容易取到1，从数值解的图形上的突变就可以看出，突变主要就是由于数值解取到1造成的。上述x用matlab计算数值解，其实取到1的极值点只有两个，对于其他的极值点我在程序中也处理了，令极大值为1，极小值为-1（任达千论文可证，上述程序这句调试代码没写），所以对于符号图形没有出现间断点甚是疑惑。（又忘了，帮我删重复贴）

yghit08

shuihai707 发表于 2013-4-23 17:16

                               
登录/注册后可看大图

数值解很容易取到1，从数值解的图形上的突变就可以看出，突变主要就是由于数值解取到1造成的。上述x用mat ...

数值解出现间断主要是在当x=1时，各种舍入误差造成无法确定是pi还是-pi，那么这样就很容易出现相位差在这一点不稳定，或大或小，造成计算瞬时频率时出现这个问题。而在做符号计算的时候，没有这些舍入误差，不会出现相位差畸变的情况。即在x=1时，那么其diff(x，'t')部分就是等于0.或者是你考察一下公式，在当x=1时是怎么处理的。
没有说明白的话，笼统的解释而且是合理的就是：一个是离散的数值解，一个是连续的理论解，只是将理论解离散了而已。刚想到的：离散点的微分表达式是基于差分方式的，即x对t的导是前后两个数据之间的差做出来的，是一种逼近，这也是说明了这种做法本身就是误差的来源，造成在x=1时，其对t的导不是0，或者小于零或者大于0的一个很小的数，那么就不是0/0型。而符号求解就不存在这个问题，在x=1时，其对t的导就是0.

shuihai707

yghit08 发表于 2013-4-23 17:32

                               
登录/注册后可看大图

数值解出现间断主要是在当x=1时，各种舍入误差造成无法确定是pi还是-pi，那么这样就很容易出现相位差在这 ...

讲的很透彻！十分感谢，得出一结论：求瞬时频率不能考虑符号函数编程了。

yghit08

shuihai707 发表于 2013-4-23 18:03

                               
登录/注册后可看大图

讲的很透彻！十分感谢，得出一结论：求瞬时频率不能考虑符号函数编程了。

如果你能将其做成模型的话，还是可以的。但是，如果你能做成模型，已经不需要符号计算来求瞬时频率了。

qqcd

感觉是acos(x)原因？可能需要考虑2pi修正吧？因为acos(x)最后结果都是[0 pi]，但是相位确实递增的，所以如果直接abs(acos(x))这样会在0和pi附近出现跳跃？

yghit08

qqcd 发表于 2013-4-24 11:46

                               
登录/注册后可看大图

感觉是acos(x)原因？可能需要考虑2pi修正吧？因为acos(x)最后结果都是[0 pi]，但是相位确实递增的，所以如果 ...

你感觉错了！这种连续变化的时候，相位是不会出现跳变的，因为x的值不会出现跳变，继而相位差不会出现跳变！

shuihai707

yghit08 发表于 2013-4-23 18:08

                               
登录/注册后可看大图

如果你能将其做成模型的话，还是可以的。但是，如果你能做成模型，已经不需要符号计算来求瞬时频率了。

我看了一些三点、五点求数值微分的程序，还是用符号计算的，这样用这些程序计算瞬时频率还是不靠谱吧，不知版主有没有比diff函数精度高的差分方法？