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

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

正常数值计算瞬时频率
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();
符号计算瞬时频率
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();

高手指点迷津!!!

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

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

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

yghit08 发表于 2013-4-23 15:39 static/image/common/back.gif
一个是连续的方法,一个是离散的方法

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

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

shuihai707 发表于 2013-4-23 16:01 static/image/common/back.gif
符号法给定的时间值是离散的啊?对于y的表达式在极值点处应该出现间断点啊?能详细说明吗?

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

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

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

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

本帖最后由 shuihai707 于 2013-4-23 16:40 编辑

yghit08 发表于 2013-4-23 16:11 static/image/common/back.gif
自己考虑一下。符号计算中,t不给离散数据你怎么在计算机中表达出来?利用符号计算的过程是连续的,亦即它 ...
原来这样,但对于公式,图形应该在极值点处出现间断点的(帮我把重复的帖子删了,操作失误,自己删不了)

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

shuihai707 发表于 2013-4-23 16:39 static/image/common/back.gif
原来这样,但对于公式,图形应该在极值点处出现间断点的(帮我把重复的帖子删了,操作失误,自己删不了) ...

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

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

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

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

yghit08 发表于 2013-4-23 17:00 static/image/common/back.gif
至于间断点,你自己分析整个过程中x能否取到1.这和你数值解的情况一样,如果数值解没有出现,那么你怎么能 ...
数值解很容易取到1,从数值解的图形上的突变就可以看出,突变主要就是由于数值解取到1造成的。上述x用matlab计算数值解,其实取到1的极值点只有两个,对于其他的极值点我在程序中也处理了,令极大值为1,极小值为-1(任达千论文可证,上述程序这句调试代码没写),所以对于符号图形没有出现间断点甚是疑惑。(又忘了,帮我删重复贴)

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

本帖最后由 yghit08 于 2013-4-23 17:38 编辑

shuihai707 发表于 2013-4-23 17:16 static/image/common/back.gif
数值解很容易取到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 发表于 2013-4-23 18:03

yghit08 发表于 2013-4-23 17:32 static/image/common/back.gif
数值解出现间断主要是在当x=1时,各种舍入误差造成无法确定是pi还是-pi,那么这样就很容易出现相位差在这 ...

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

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

shuihai707 发表于 2013-4-23 18:03 static/image/common/back.gif
讲的很透彻!十分感谢,得出一结论:求瞬时频率不能考虑符号函数编程了。

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

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

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

yghit08 发表于 2013-4-24 11:55

qqcd 发表于 2013-4-24 11:46 static/image/common/back.gif
感觉是acos(x)原因?可能需要考虑2pi修正吧?因为acos(x)最后结果都是,但是相位确实递增的,所以如果 ...

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

shuihai707 发表于 2013-4-24 12:40

yghit08 发表于 2013-4-23 18:08 static/image/common/back.gif
如果你能将其做成模型的话,还是可以的。但是,如果你能做成模型,已经不需要符号计算来求瞬时频率了。

我看了一些三点、五点求数值微分的程序,还是用符号计算的,这样用这些程序计算瞬时频率还是不靠谱吧,不知版主有没有比diff函数精度高的差分方法?
页: [1] 2
查看完整版本: 符号和数值计算瞬时频率的差别