24小时热门版块排行榜    

查看: 1965  |  回复: 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 的主题更新
☆ 无星级 ★ 一星级 ★★★ 三星级 ★★★★★ 五星级
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] filecode +6 documentary 2026-08-10 6/300 2026-08-10 17:09 by wdc28
[基金申请] 我的国基提前知道中了,可是同事的操作让我实在接受不了,怎么会有这样的人 +8 家与远方 2026-08-10 11/550 2026-08-10 17:09 by 六两废铜
[基金申请] FileCode能看出啥? +4 要乐观耀哥 2026-08-10 8/400 2026-08-10 15:02 by 天神眷顾
[基金申请] 应该是93bebmhtak前后十一个字符比较关键 +14 Lanmanbaby 2026-08-09 22/1100 2026-08-10 14:57 by wwncly
[基金申请] 8月时间戳变的,举个手。玩一下,释放压力 +11 archvillain 2026-08-04 13/650 2026-08-10 14:41 by seekfind
[基金申请] 好奇怪的filecode +4 布布和一二 2026-08-08 5/250 2026-08-10 14:20 by 冰心玉壶晴
[基金申请] 面上项目filecode邪修 +5 西山十月 2026-08-09 7/350 2026-08-10 07:32 by 仁砚薪传
[基金申请] fileCode有新解读? +10 Tide man 2026-08-08 18/900 2026-08-09 12:55 by 仁砚薪传
[基金申请] 关于filecode,很负责任的告诉大家 +6 爱看书的可乐 2026-08-08 7/350 2026-08-08 22:13 by a_niu
[基金申请] 关于豆爷回答的JTJC与%2F数量 +5 yang182083 2026-08-06 7/350 2026-08-08 18:29 by zhanghaozhu
[基金申请] 关于filecode +4 布布和一二 2026-08-07 7/350 2026-08-07 22:55 by zhanghaozhu
[基金申请] 娱乐 +6 Tide man 2026-08-03 6/300 2026-08-07 22:40 by 铁帽子农民
[基金申请] 固定端突然变了,今天 +6 archvillain 2026-08-06 10/500 2026-08-07 16:03 by 医学老男孩
[基金申请] 听说今天filecode变了 +24 布布和一二 2026-08-06 47/2350 2026-08-07 16:02 by zhiyanjiang
[基金申请] filecode变化情况 +6 布布和一二 2026-08-07 22/1100 2026-08-07 14:45 by 且听虎啸
[基金申请] 大家散了吧,后缀研究没有意义,别浪费时间了,过好目前的每一天,不要焦虑 +5 Tide man 2026-08-06 7/350 2026-08-07 13:11 by 医学老男孩
[基金申请] filecode +8 布布和一二 2026-08-06 11/550 2026-08-06 20:41 by tangpu318
[基金申请] 求各位大神看下 100+6 hpkpkpkp 2026-08-05 33/1650 2026-08-06 14:49 by zhiyanjiang
[基金申请] 影响面上的因素 +8 布布和一二 2026-08-05 11/550 2026-08-06 10:41 by 宝贝虫子
[基金申请] 纯娱乐,不喜欢勿喷 +7 Tide man 2026-08-04 10/500 2026-08-04 15:10 by loufangrui
信息提示
请填处理意见