声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2254|回复: 6

[综合讨论] 请教: 样条插值积分

[复制链接]
发表于 2006-8-14 11:17 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
我想对一列有误差的加速度值进行样条函数插值积分以求得位移(使用已知解析表达式的函数是为了验证用此方法的正确性),以下为我的程序,但是结果不正确,请各位高手帮忙看看,我的思路有问题还是函数用的不对?多谢多谢!
n=1024;  
t=0:1/fs:n/fs
y1=50*sin(4*pi*t)+100*sin(10*pi*t);    %各采样得到的加速度值
sp1=csapi(t,y2);                                 %三次样条插值
Cformint=fnint(sp1,0);       %对插值得到的函数进行积分
varCformint=fnval(Cformint,t);       %得到积分后在各时间点的值(速度值)
sp2=csapi(t,varCformint);      %对插值得到的速度值再次进行三次样条插值
Cformint2=fnint(sp2,0);      %对插值得到函数进行积分以求得位移
varCformint2=fnval(Cformint2,t);     %得到积分后位移在各时间点的值
figure
plot(t,varCformint2)

得到的结果不是正弦曲线而是一条斜率为正的直线,这是为什么呢,请各位高手帮忙解答一下,感激不尽!
回复
分享到:

使用道具 举报

发表于 2006-8-14 21:42 | 显示全部楼层
估计是积分造成的,希望给出完整的代码
 楼主| 发表于 2006-8-15 15:12 | 显示全部楼层
实在不好意思,我为了尽量减少帖子的字数,删减了原来程序的一部分内容,没想到把有用的部分也给删掉了,以下的程序是可以运行的:

n=1024;  
fs=20;
t=0:1/fs:n/fs;
y1=100*sin(10*pi*t);   
sp1=csapi(t,y1);                                 
Cformint=fnint(sp1,0);      
varCformint=fnval(Cformint,t);      
sp2=csapi(t,varCformint);      
Cformint2=fnint(sp2,0);     
varCformint2=fnval(Cformint2,t);   
figure
plot(t,varCformint2)

我的问题是y1这个正弦函数如果进行两次积分的话,得到的结果应该还是正弦函数,
但是我对其用样条函数进行插值,然后进行数值积分,得到却不是正弦函数曲线,自己找不出原因,请各位高手多加指点,多谢多谢!
发表于 2006-8-15 17:19 | 显示全部楼层
Cformint=fnint(sp1,0);

是这个语句造成的
实际积分的结果应该是F(x)+C,你用这个语句表示的是C=0
而实际上C并不等于零,换句话说就是在你第一次积分的过程中引入了直流量

所以第二次积分的时候这个直流量就是就变成了斜线量

由于这个量在数值上远大于波动量,所以看上去是斜线,实际上放大看的话也应该是波动的

这就需要将程序中将着一直流量滤掉
 楼主| 发表于 2006-8-17 09:53 | 显示全部楼层
非常感谢您的解答!  能否再请问您该如何去除这个直流项?
因为我既使将   Cformint=fnint(sp1,0);   这个语句改成    Cformint=fnint(sp1,-100/10/pi); 其中,-100/10/pi是当t=0时使位移为0,得到的积分常数C的值, 再运行程序,得到的结果仍然是斜直线。 这是怎么回事呢?请多多指教!
发表于 2006-8-17 17:07 | 显示全部楼层
一般采用滤波的方法,这个问题论坛讨论过好多

比如:http://forum.vibunion.com/forum/viewthread.php?tid=20424
http://forum.vibunion.com/forum/viewthread.php?tid=329

[ 本帖最后由 happy 于 2006-8-17 17:08 编辑 ]
 楼主| 发表于 2006-8-18 10:07 | 显示全部楼层

谢谢!

谢谢您的回复,我会找相关讨论来好好学习的!再次感谢!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-5-19 00:14 , Processed in 0.050270 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表