打印

要用DSP做一个频响函数,IIR系统的频率响应算法可不可以用

[复制链接]
665|1
手机看帖
扫描二维码
随时随地手机跟帖
跳转到指定楼层
楼主
qaz098xsw|  楼主 | 2016-6-7 22:23 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
我做的项目是用激振锤锤击悬臂梁,DSP采集激振锤返回的力信号的电压值和安装在悬臂梁上的加速度传感器返回的加速度信号的电压值,将这两路信号做频响函数,但不知道这个程序能不能应用,各位大神有空的帮忙看看哈,万分感谢。我采样频率设置的是6K,采样点数是512个点。
程序如下:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#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<n;k++)           
         {
                 sf=s*k;            
                 bsum.real=b[0];bsum.imag=0.0;            
                 if(lb>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;         
         }

相关帖子

沙发
Sode| | 2016-6-19 15:41 | 只看该作者
支持一下

使用特权

评论回复
发新帖 我要提问
您需要登录后才可以回帖 登录 | 注册

本版积分规则

632

主题

842

帖子

3

粉丝