微波EDA网,见证研发工程师的成长!
首页 > 研发问答 > 嵌入式设计讨论 > DSP学习交流 > 要用DSP做一个频响函数,网上找到了一个IIR系统的频率响应算法实例,有点看不太懂我是不是能使用

要用DSP做一个频响函数,网上找到了一个IIR系统的频率响应算法实例,有点看不太懂我是不是能使用

时间:10-02 整理:3721RD 点击:
我做的项目是用激振锤锤击悬臂梁,DSP采集激振锤返回的力信号的电压值和安装在悬臂梁上的加速度传感器返回的加速度信号的电压值,将这两路信号做频响函数,但不知道这个程序能不能应用,各位大神有空的帮忙看看哈,万分感谢。我采样频率设置的是6K,采样点数是512个点。
程序如下:
#include
#include
#include
#include "msp.h"
void  miirres(float a[],float b[],int lb,int la,complex h[],int n)
{
/*---------------------------------------------------------------------
Routinue miirres:to obtain the frequence responce (Amplitude and phase)     
     H(exp(jw)) of system  function H(z)=B(z)/A(z);            
                 b(0)+b(1)z^(-1)+ ... + b(lb)z^(-lb)      
        H(z)=------------------------------------            
                 1.0 +a(1)z^(-1)+ ... + a(la)z^(-la)
input parameters:
a :la+1 dimensioned real array. a(0)=1.0,   
   a(n) (n=1 to la) is the coef. of z^(-n).
b :lb+1 dimensioned real array. b(n) (n=0 to la) is the coef. of z^(-n).
n :data point number of frequency response from 0.0 to pi.the frequency
   resolution is pi/n.
output parameters:
H :n dimensioned complex array. frequency response is stored in   
   H(0) to H(n-1), which is normalized from 0.0 to pi.
Note: If H(k)=0/0,then let H(K)=H(k+1) (This is for the case of that      
      H(exp(jw)) is a sinc function)     
: If la=0: for FIR system.         
                                in Chapter 2
---------------------------------------------------------------------*/        
        complex z,bsum,asum;        
        int k,i;        
        int isign;        
        float pole,peak,pi,s,sf,sa,sb,temp;  
        pole=1.e-10;        
        peak=1.e+10;        
        pi=4.*atan(1.);        
        s=pi/(float)(n);        
        isign=0;        
        for(k=0;k 0)               
                        for(i=1;i<=lb;i++)                  
                                {                  
                                sb=sf*i;                  
                                z.real=cos(sb);z.imag=-sin(sb);
                                bsum.real+=b*z.real;         
                                bsum.imag+=b*z.imag;                    
                                }            
                asum.real=1.;            
                asum.imag=0.0;            
                if(la!=0)               
                        for(i=1;i<=la;i++)                  
                                {
                                sa=sf*i;                  
                                z.real=cos(sa);z.imag=-sin(sa);                  
                                asum.real+=a*z.real;                  
                                asum.imag+=a*z.imag;                    
                                }            
                if(mabs(asum)<pole)              
                        {               
                        if(mabs(bsum)<pole)              
                        isign=1;              
                        else              
                        {               
                        h[k].real=peak;              
                        h[k].imag=0.0;               
                        }               
                        }            
                        else              
                        {               
                        temp=pow(mabs(asum),2.0);              
                        h[k].real=(bsum.real*asum.real+bsum.imag*asum.imag)/temp;              
                        h[k].imag=(bsum.imag*asum.real-bsum.real*asum.imag)/temp;              
                        if(isign==1)              
                        {               
                        h[k-1]=h[k];              
                        isign=0;               
                        }               
                        }         
                }        
        return;         
        }


已退回10积分

Copyright © 2017-2020 微波EDA网 版权所有

网站地图

Top