24小时热门版块排行榜    

查看: 2004  |  回复: 0

baobiao007

木虫 (职业作家)


[资源] 【分享】最小平方反滤波(反褶积)【无重复】

//最小平方反滤波,学习反褶积必须掌握的第一个程序

#include
#include
#include
/*解Teoplize矩阵的莱文森递推算法
  t[]--矩阵的第一行数据,n个
  b[]--线性方程组最右边的一列数据
  x[]--计算结果
*/
int atlvs(double t[],int n,double b[],double x[])
  { int i,j,k;
    double a,beta,q,c,h,*y,*s;
    s=(double *)calloc(n,sizeof(double));
    y=(double *)calloc(n,sizeof(double));
    a=t[0];
    if (fabs(a)+1.0==1.0)
      { free(s); free(y); printf("fail\n"; return(-1);}
    y[0]=1.0; x[0]=b[0]/a;
    for (k=1; k<=n-1; k++)
      { beta=0.0; q=0.0;
        for (j=0; j<=k-1; j++)
          { beta=beta+y[j]*t[j+1];
            q=q+x[j]*t[k-j];
          }
        if (fabs(a)+1.0==1.0)
          { free(s); free(y); printf("fail\n"; return(-1);}
        c=-beta/a; s[0]=c*y[k-1]; y[k]=y[k-1];
        if (k!=1)
          for (i=1; i<=k-1; i++)
            s=y[i-1]+c*y[k-i-1];
        a=a+c*beta;
        if (fabs(a)+1.0==1.0)
          { free(s); free(y); printf("fail\n"; return(-1);}
        h=(b[k]-q)/a;
        for (i=0; i<=k-1; i++)
          { x=x+h*s; y=s;}
        x[k]=h*y[k];
      }
    free(s); free(y);
    return(1);
  }

//互相关
void corre(double x[],int m,double y[],int n,double z[],int l)
{
        int i,j,k;
        for(i=0; i         {
                z=0.0;
                k=i-(n-1);
                for(j=0; j                         if(j-k>=0&&j-k                                 z+=x[j]*y[j-k];
        }
}
/*  线性卷积
    y(n)=x(n)*h(n)
    m--length of x(n);
        n--length of h(n);
        l=m+n-1  length of y(n)
*/
void conv(double x[],int m,double h[],int n,double y[],int l)
{
        int i,j,k;
        for(i=0; i         {
                y=0.0;
                for(j=0; j                 {   k=i-j;
                        if(k>=0 && k                                 y += x[j]*h[k];
                }
        }
}
//产生零相位子波,num--length,f0--frequency
void WaveLet(double x[],int num,double f0)
{
        double pi=3.1415926;
        double det=0.002;//2ms采样
        double m=2.0;//余弦波衰减比例
        double a;
        a=2*f0*f0*log(m);
        for(int i=0; i                 x=exp(-a*i*i*det*det)*cos(2*pi*f0*i*det);
}

void main()
{
        const int N=129;
        const int M=2*N-1;
        double x[N];//子波
        double y[N];//计算出的反子波
        double dirac[N]={1.0};//单位冲激序列
         double r1[M], r2[N]={0.0};
        int i;

        WaveLet(x,N,40.0);
        //子波自相关
        corre(x,N,x,N,r1,M);
        for(i=0; i                 r2=r1[M/2+i];
        //形成并解Toeplize方程
         atlvs(r2,N,dirac,y);
        //反子波与子波卷积变回单位冲激序列,说明计算正确
        conv(x,N,y,N,r1,M);
        for(i=0; i             printf("r1[%d]=%f\n",i,r1);
}
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

智能机器人

Robot (super robot)

我们都爱小木虫

找到一些相关的精华帖子,希望有用哦~

科研从小木虫开始,人人为我,我为人人
相关版块跳转 我要订阅楼主 baobiao007 的主题更新
☆ 无星级 ★ 一星级 ★★★ 三星级 ★★★★★ 五星级
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +4 FChKX8meu1ds 2026-11-09 8/400 2026-11-14 22:04 by Equinoxhua
[考博] 27博士 +6 锦鲤送福至 2026-11-10 6/300 2026-11-14 21:23 by xupeii
[公派出国] 售SCI一区文章,我:8O5.5.1.O5.4,科目全,可伽急 +3 TZqyNTBa07NP 2026-11-13 3/150 2026-11-14 16:09 by KXypaGgBBlw7
[硕博家园] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +6 zSbvVX9JFl9J 2026-11-09 9/450 2026-11-13 11:13 by Equinoxhua
[考研] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +8 zSbvVX9JFl9J 2026-11-09 9/450 2026-11-13 11:11 by Equinoxhua
[硕博家园] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +5 C6CRn66RYYXS 2026-11-10 9/450 2026-11-13 11:02 by Equinoxhua
[考博] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +4 cdqGD6ftHmbl 2026-11-11 5/250 2026-11-12 12:28 by ECLTXnJnRSXf
[教师之家] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +4 Z7qfCsdIiD1T 2026-11-10 7/350 2026-11-12 11:22 by HajghKOMynmz
[硕博家园] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +5 Z7qfCsdIiD1T 2026-11-10 8/400 2026-11-12 11:20 by HajghKOMynmz
[考博] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +4 czAiolzR0IJC 2026-11-10 7/350 2026-11-12 11:11 by HajghKOMynmz
[教师之家] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +5 czAiolzR0IJC 2026-11-10 7/350 2026-11-12 11:09 by HajghKOMynmz
[硕博家园] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +4 czAiolzR0IJC 2026-11-10 8/400 2026-11-12 10:49 by HajghKOMynmz
[公派出国] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +7 FChKX8meu1ds 2026-11-09 9/450 2026-11-12 10:03 by HajghKOMynmz
[考博] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +5 FChKX8meu1ds 2026-11-09 8/400 2026-11-12 09:46 by HajghKOMynmz
[教师之家] 售SCI一区T0P文章,我:8.O.55.1.O.54,科目齐全,可+急 +3 TgvkbZZl9tqB 2026-11-11 4/200 2026-11-12 08:15 by HajghKOMynmz
[考研] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 cdqGD6ftHmbl 2026-11-11 3/150 2026-11-12 07:53 by HajghKOMynmz
[教师之家] 售SCI文章,我:8O5.5.1.O.54,科目齐全,可+急 +4 C6CRn66RYYXS 2026-11-10 4/200 2026-11-12 06:03 by 9fOTqj074tIg
[基金申请] 售SCI一区T0P文章,我:8.O.55.1.O.54,科目齐全,可+急 +5 FChKX8meu1ds 2026-11-09 7/350 2026-11-12 05:30 by 9fOTqj074tIg
[教师之家] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +4 FChKX8meu1ds 2026-11-09 7/350 2026-11-12 05:19 by 9fOTqj074tIg
[基金申请] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 FChKX8meu1ds 2026-11-09 7/350 2026-11-12 05:19 by 9fOTqj074tIg
信息提示
请填处理意见