基于DSP的语音处理系统的设计

摘要近年来，随着DSP技术的普及和低价格、高性能DSP芯片的出现，DSP已越来越多地被广大的工程师所接受越来越广泛地被应用于各个领域，并且已日益显示出其巨大的优越性。DSP是利用专门或通用的数字信号处理芯片，以数字计算的方法对信号进行处理，具有处理速度快、灵活、精确、抗干扰能力强、体积小及可靠性高等优点，满足了对信号快速、精确、实时处理及控制的要求。本次设计基于TLV320AIC23和TMS320VC5416两种芯片设计并实现了一种语音录音、语音编码、语音解码、语音处理和回放的系统。通过软件和硬件结合对该系统进行设计，使本次设计的语音处理系统具有强大的数据处理能力并配有灵活的接口电路，可以作为一种语音信号处理算法研究和实时实现的通用平台，对语音编码在DSP上的实时实现进行了简单的研究，从而掌握了算法移植的一般流程，为能够在高速DSP硬件平台设计及系统应用开发方面取得成功奠定基础。　　关键词：DSP；数据采集； TLV320AIC23；TMS320VC5416。目录摘要 I 第1章绪论 1 1.1 DSP的发展及应用 1 1.2 语音信号处理系统概述 2 第2章 DSP芯片介绍 3 2.1 TLV320AIC23简介 3 2.2 TMS320VC5416简介 3 第3章系统设计 4 3.1系统硬件设计 4 3.1.1系统结构框图 4 3.1.2 DSP处理器 5 3.1.3 A/D电路 5 3.1.4 D/A电路 7 3.2系统软件设计 10 3.2.1 TMS320VC5416初始化 10 3.2.2 TLV320AIC23初始化 10 第4章总结 11 参考文献 12 致谢 13 附录 14 第1章绪论近年来，在数字信号处理领域有着绝对优势的DSP技术得到了迅速发展，不仅在通信计算机领域大显身手，并已逐渐渗透到人们日常消费领域。正因为如此，DSP应用越来越得到普遍重视。 DSP作为可编程数字信号处理专用芯片是微型计算机发展的一个重要分支，也是数字信号处理理论实用化过程的重要技术工具。DSP器件分为两大类：一类是专门用于FFT、FIR滤波、卷积等运算的芯片，称为专用DSP器件；另一类是可以通过编程完成各种用户要求的信息处理任务的芯片，称为通用数字信号处理器件。 1.1 DSP的发展及应用最初的DSP器件只是被设计用以完成复杂数字信号处理算法。这可以追溯到20世纪50年代到60年代，那时数字信号处理技术刚刚起步。由于一般的数字信号处理算法运算量大，因此，算法只能在大型计算机上进行模拟仿真，无法实现数字信号处理。60年代中期，快速傅里叶算法的出现及大规模集成电路的发展，奠定了硬件完成数字信号处理算法和数字信号处理理论实用化的重要技术基础，从而促进了近40年来DSP技术与器件的飞速发展。通用DSP器件的发展可分为三个阶段：第一阶段（1980年前后），DSP雏形阶段。第二阶段（1990年前后），DSP的成熟阶段。第三阶段（2000年以后），DSP的完善阶段。目前，DSP的发展非常迅速。硬件结构方面主要是向多处理器的并行处理结构、便于外部数据交换的串行总线传输、大容量片上RAM和ROM、程序加密、增加I/O驱动能力、外围电路内装化、低功耗等方面发展。软件方面主要是综合开发平台的完善，使DSP的应用开发更加灵活方便。目前，DSP芯片的价格越来越低，性能价格比日益提高，具有巨大的应用潜力。DSP芯片的主要应用： ① 信号处理——数字滤波，自适应滤波，快速傅里叶变换，相关运算，频谱分析，卷积，波形产生等； ② 语音处理——语音编码，语音识别，语音合成，文本—语音转换等； ③ 图象图形处理——三维图形转换，机器人视觉，图象转换及压缩，模式识别，图象增强等； ④ 控制——司服控制，机器人控制，自适应控制，神经网络控制等； ⑤ 军事——保密通信，雷达及声音信号处理，导航及制导，调制解调，全球定位，搜索与跟踪等； ⑥ 仪器仪表——频谱分析，函数发生器，模态分析，暂态分析等； ⑦ 通讯——回音相消，高速调制解调器，数字编码与解码，自适应均衡，移动电话，扩展通讯，噪音对消，网络通讯等； ⑧ 消费电子——高清晰度电视，音乐合成器，智能玩具，游戏等； ⑨ 医学——助听器，病员监控，超声波设备，自动诊断设备，胎儿监控等。 1.2 语音信号处理系统概述语音处理在现代通信中应用非常广泛，主要有语音编码、语音识别、语音合成、语音邮件、语音存储等。典型的语音处理系统如下图1.1所示：图1.1 典型的语音处理系统图中的输入信号可以有各种各样的形式。例如，它可以是麦克风输入的语音信号或是电话线已调的数据信号，可以是编码后在数字链路上传输或存储在计算机里德摄像机图象信号等。输入信号首先进行带限滤波和抽样，然后进行A/D变换将信号变成数字比特流。根据奈奎斯特抽样定理，为保证信息不丢失，抽样频率至少是输入带限信号最高频率的2倍。 DSP芯片的输入是A/D变换后得到的以抽样形式表示的数字信号。DSP芯片对输入的数字信号进行某种处理。数字处理是DSP系统的关键，这与其它系统（如电话交换系统）有很大的不同。在交换系统中，处理器的作用是进行路由选择，它并不对输入数据进行修改。因此两者虽然都是实时系统，但两者的实时约束条件却有很大不同。最后，经过处理后的数字样值再经D/A变换转换为模拟样值。之后进行内插和平滑滤波就会得到连续的模拟波形。上面给出的典型的DSP语音处理系统，根据不同的用途应有不同的变动。第2章 DSP芯片介绍 2.1 TLV320AIC23简介 TLV320AIC23（简称AIC23）是TI公司的一款高性能Codec芯片。主要特性有：内置耳机输出放大器，支持MIC和LINE IN两种输入方式（二选一）。且对输入和输出都具有可编程增益调节；芯片中的A/D转换器和D/A转换器采用多位的Sigma-Delta技术，数据传输字长为16、20、24、32bit，采样率为8kHz ~96kHz；在采样率为96kHz情况下A/D转换器信噪比达到90dB，D/A转换器达到100dB；回放模式下功率为23mW，省电模式下更是小于15uW；只占用25mm的面积。基于上述优点，AIC23是可移动的数字音频播放和录音使用中的模拟输入输出等应用系统的理想选择，例如MP3播放器等。 2.2 TMS320VC5416简介 TMS320VC5416（以下简称VC5416）是TI公司的一款16bit定点高性能DSP，是TMS320VC54x系列中的第3代芯片。主要特性有：速率最高达160MI/s；3条16bit数据存储器总线和1条程序存储器总线；1个40bit桶形移位器和2个40bit累加器；1个17×17乘法器和1个40bit专用加法器；最大8M×16bit的扩展寻址空间，内置128k×16bit的RAM和16k×16bit的ROM；3个多通道缓冲串口（McBSP）；配有PCM3002，可对语音进行A/D和D/A转换。由于VC5416功耗低，性能高，其分开的数据和指令空间使该芯片具有高度的并行操作能力，在单周期内允许指令和数据同时存取，再加上高度优化的指令集，使得该芯片具有很高的运算速度并且该芯片本身具有丰富的片内存储器资源和多种片上外设，因此在工程界得到广泛应用，尤其是在语音编码和通信应用方面。第3章系统设计 3.1系统硬件设计 3.1.1系统结构框图音频系统应该具有较宽的动态范围，选择16~24位的ADC和DAC能完全捕获或恢复高保真的音频信号。系统的核心芯片(DSP)选用美国TI公司的TMS320VC5402[1](以下简称C5402)。 DSP芯片模块是整个实时语音处理系统的核心部分，它对经数字化的信号进行压缩，编解码等。 A/D转换模块功能是把模拟信号数字化，包括采集和量化，这部分为DSP处理语音数字信号做好了准备；D/A转换模块就是把数字信号转换为模拟的信号，输出音频信号。 SDRAM（动态随机存储器）存储器模块主要是为DSP处理器扩展存储容量，达到要求的存储容量；但要注意的是要与DSP处理器的速度相匹配，以便良好的运行。电源模块是为内部芯片及周边系统电路提供能量的部分。系统结构框图图3.1如下所示：图3.1 系统结构框图 3.1.2 DSP处理器作为DSP家族高性价比代表的16位定点DSP芯片，C5402适用于语音通信等实时嵌入应用场合。与其它C54X芯片一样，C5402具有高度灵活的可操作性和高速的处理能力。其性能特点如下：操作速率可达100MIPS；具有先进的多总线结构，三条16位数据存储器总线和一条程序存储器总线；40位算术逻辑单元(ALU)，包括一个40位桶形移位器和两个40位累加器；一个17×17乘法器和一个40位专用加法器，允许16位带/不带符号的乘法；整合维特比加速器，用于提高维特比编译码的速度；单周期正规化及指数译码；8个辅助寄存器及一个软件栈，允许使用业界最先进的定点DSP C语言编译器；数据/程序寻址空间为1M×16bit，内置4K×16bit ROM和16k×16bit RAM；内置可编程等待状态发生器、锁相环(PLL)时钟产生器、两个多通道缓冲串口、一个与外部处理器通信的8位并行HPI口、两个16位定时器以及6通道DMA控制器且低功耗。与C54X系列的其它芯片相比，5402具有高性能、低功耗和低价格等特点。它采用6级流水线，且当RPT(重复指令)时，一些多周期的指令就变成了单周期的指令；芯片内部RAM和ROM可根据PMST寄存器中的OVLY和DROM位灵活设置。这些都有利于算法的优化。 C5402采用3.3V和1.8V电源供电，其中I/O采用3.3V电源供电，芯片的核采用1.8V电源供电。而实际常用的只有5V电源，所以必须采用电源转换芯片。选用TPS7301和TPS7333两块电源转换芯片(它们都是TI公司为配合DSP而设计的电源转换芯片)，分别接上适当的外围电阻，构成电阻分压器，即可调整两块芯片的输出电压分别为3.3V和1.8V。 3.1.3 A/D电路 PCM1800是双声道单片Δ－Σ型20位ADC，单+5V电源供电，信噪比为95dB，动态范围为95dB，其内部嵌有高通滤波器，具有PCM音频接口和四种数据格式，分为主控和受控两种模式，采样频率可选为32kHz、44.1KHz和48KHz。 PCM1800构成音频信号采集系统时，主要涉及到BCK(位时钟信号)、LRCK(采样时钟信号)、FSYNC(帧同步信号)、DOUT(数字信号输出)、SYSCLK(系统时钟输入)这几个对时序有要求的引脚。通过对引脚MODE0和MODE1进行编程，可让PCM1800工作于主控模式(Master Mode)。此时，BCK、LRCK、FSYNC均作为输出，其时序由PCM1800内部的时钟产生电路控制。但SYSCLK只能由外部提供(这里用C5402的TOUT脚输出信号提供)。 PCM1800的系统时钟只能是256fs、384fs或者512fs，这里fs是音频信号采样频率。在主控模式时，FSYNC用来指明PCM1800的DOUT输出的有效数据，它的上升沿表明一帧数据的起始，下降沿表明一帧数据的结束。FSYNC的频率是采样时钟频率LRCK的2倍。在此模式下，位时钟信号BCK的频率是采样时钟频率LRCK的64倍。通过对PCM1800的FMT0、FMT1两引脚编程(FMT0＝1，FMT1＝0)，可以设置PCM1800输出的数据格式为20位的IIS格式。为了保证在数据处理时不影响新数据的接收以及在接收数据时不中断正在进行的数据处理过程，采用了多通道缓冲同步串口(McBSP)。 PCM1800与C5402连接后，C5402使用缓冲串口0接收数据，各种同步信号由PCM1800产生，C5402是被动接收各种信息。PCM1800与C5402的硬件接线图如图3.2所示。图3.2 PCM1800与C5402的硬件接线图电源管理功能模块所用器件: TPS73HD301( 3.3-V/Adjustable Output 该芯片一端输入可调,范围是（1.2-9.7V）电源模块管脚图3.3如下所示：图3.3 电源模块管脚图 3.1.4 D/A电路 PCM1744是双声道立体声DAC，包含数字滤波器和输出放大器，动态范围为95dB，具有多种采样频率可选，最高可达96kHz。采用24位的IIS数据输入格式。PCM1744的操作主要涉及到LRCIN(采样时钟信号输入)、BCKIN(位时钟信号输入)、SCKI(系统时钟输入)、DIN(数据输入)这几个对时序有要求的引脚。PCM1744与C5402连接后，C5402使用缓冲串口1发送数据，各种时钟信号均由C5402产生，PCM1744被动接收各种信息。PCM1744的系统时钟信号(SCKI)由C5402的TOUT引脚提供，TOUT是C5402的定时器输出信号引脚，有较强的驱动能力，可以驱动多个芯片。PCM1744的数据接收时钟格式必须是IIS格式，C5402在缓冲串口寄存器中设置各种时钟方式时，必须满足IIS格式的要求。C5402作为主动工作器件，可以对其缓冲串口输出信号进行调整。输出的采样时钟信号、位时钟信号可以在McBSP寄存器SRGR1和SRGR2中设置，设置遵循图3.4的原则。图3.4 C5402时钟发生流程图基本的时钟信号可以来自CPU时钟，也可以来自晶振时钟，这在SRGR2寄存器中的第13位设置。基本时钟输入后，经CLKGDV(SRGR1的第7位到第0位)所设置的值进行第一次分频，得到位时钟信号(由BCLKX1脚输出)。值得注意的是，位时钟信号最高为DSP频率的一半。位时钟信号经FPER(SRGR2的第11位到第0位)和FWID(SRGR1的第15位到第8位)所设置的值进一步分频得到采样时钟信号(由BFSX1脚输出)，FPER和FWID分别设置采样时钟信号的低电平和高电平的时间值。C5402与PCM1744的硬件接线如图3.5所示。图3.5 PCM1744与C5402接线图 PCM1800完成音频信号采集后，在DSP的外扩程序存储器中嵌入相应的处理算法，语音信号经处理后，再从PCM1744输出。复位电路:所用芯片为74HC14。复位电路图3.6如下所示：图3.6 复位电路图存储器模块所用芯片为：MT48LC8M8A2TG-75、存储容量Density 为64Mb、数据宽度16位、工作电压 3.3V、 TSOP封装 54管脚、时钟速率133 MHz、存储器模块图3.7如下所示：图3.7 存储器模块图 3.2系统软件设计 3.1.1 TMS320VC5416初始化　　 SWWSR=0x7fff;/程序、数据、I/O空间　　SWCR=0x0001;/等待周期为7×2=14 　　BSCR=0x8006;/按32KW分区，HD[7:0]，D[15:0]Hold 　　CLKMD=PLL_DIV_INIT; 　　Waitloop(0x0400); 　　CLKMD=PLL_LOCK_INIT_X(5); 　　Waitloop(0x0400); 　　PMST=0x0168;/中断向量表定位在0X100，MP/MC=1 　　OVLY=1 DROM=1; 3.2 TLV320AIC23初始化为使AIC23正常工作并产生预期的音频效果，必须对其相应的寄存器进行配置。首先对VC5416的I2C模块初始化，将AIC23总线上的地址写入从机地址寄存器ICSAR；再把相应的AIC23内部映射寄存器的地址和待写数据合并为16bit控制字，逐次写入ICDXR，并通过I2C总线发送给AIC23，即可完成对AIC23的初始化配置。　 AIC23初始化的部分源代码： Unsigned int codec_buf[9]={OX1e00，OXOc00，OxO81a OxOaO4，OxOe01，0x1020，0x1021，0x0117，OxO5f9}; 　　Port_sub_address=(unsigned int*)MCBSP_SPSA_ADDR(1); 　　Port_sub_index_reg=(unsigned int*)MCBSP_SPAD_ADDR(1); 　　*Port_sub_address=MCBSP_SPCR2_SUBADDR; 　　For(i=O;i<9;i++){ 　　While(!(*Port_sub_index_reg&(MASK_BIT(XRDY)))); 　　Set_codec_cs_low(); 　　MCBSP1_DXR1=codec_buf[i]; 　　While(!(*Port_sub_index_reg&(MASK_BIT(XRDY)))); 　　Set_codec_cs_high();} 第4章总结经过这次的课程设计，我发现我在DSP这方面学得不够，很多东西都学得不够全面，掌握得不够深，不能熟练地把它们应用在实践当中。这次在刘伟春老师的细心指导和同学的热心帮助，以及自己上网查找资料下，还算比较顺利地完成了本次课程设计的任务。这次课程设计使我对DSP方面的知识有更深的理解，强化了自己的基础知识，也深刻体会到DSP技术应用领域的广泛。同时对CCS集成环境更为熟悉了，为我在今后的工作中奠定了坚实的实践基础。通过这次课程设计让我明白基础知识的重要性，同时也要理解更多的有关它的知识，并且很好地运用到实践当中，也让我知道了要好好地学习，不能懈怠。参考文献 [1]戴明桢等编著.TMS320C54X DSP 结构原理及应用. 北京：航空航天大学出版社，第2版，2007； [2]彭启琮编著.DSP技术的发展与应用.北京：高等教育出版社，2002； [3]胡广书编著.数字信号处理理论、算法与实现.北京：清华大学出版社，2005； [4]张雄伟，曹铁勇．DSP芯片的原理与开发应用(第二版)[M]．北京：电子工业出版社，2000； [5]郝软层，徐金甫.基于DSP芯片的MELP声码器的算法实现[J].微计算机信息，2006； [6]任丽香，马淑芬，李方慧.TMS220600系列DSP的原理与应用[M]．北京：电子工业出版社，2000； [7]北京合众达电子技术有限公司编著.SEED-DTK系列实验手册.北京合众达电子技术有限公司出版，2007。致谢在这次课程设计过程中，我要感谢每一个帮助过我的人。本论文是在刘老师的悉心指导下完成的，刘老师对我的论文提出了很多宝贵的意见，帮助我开拓研究思路，精心点拨、热忱鼓励。同时，刘老师渊博的学识、严谨的治学态度也令我十分敬佩，是我以后学习和工作的榜样。在整个设计过程中我懂得了许多东西，也培养了我独立工作的能力，树立了对自己工作能力的信心，相信会对今后的学习工作生活有非常重要的影响，同时，也让我知道了那些基础知识的重要性。本论文的顺利完成，离不开我们老师悉心教导、同学和朋友的关心和帮助。总之，感谢每一位关心过我，爱护过我的人。最后，再次感谢各位老师各位同学的帮助和支持，衷心地谢谢你们！附录源代码如下： #include void iirbcf(ifilt,band,ns,n,f1,f2,f3,f4,db,b,a) double b[],a[],f1,f2,f3,f4,db; /* TMS320VC5416初始化*/ SWWSR=0x7fff;/程序、数据、I/O空间 SWCR=0x0001;/等待周期为7×2=14 BSCR=0x8006;/按32KW分区，HD[7:0]，D[15:0]Hold CLKMD=PLL_DIV_INIT; Waitloop(0x0400); CLKMD=PLL_LOCK_INIT_X(5); Waitloop(0x0400); PMST=0x0168;/中断向量表定位在0X100，MP/MC=1 OVLY=1 DROM=1; /* TLV320AIC23初始化*/ Unsigned int codec_buf[9]={OX1e00，OXOc00，OxO81a OxOaO4，OxOe01，0x1020，0x1021，0x0117，OxO5f9}; Port_sub_address=(unsigned int*)MCBSP_SPSA_ADDR(1); Port_sub_index_reg=(unsigned int*)MCBSP_SPAD_ADDR(1); *Port_sub_address=MCBSP_SPCR2_SUBADDR; For(i=O;i<9;i++){ While(!(*Port_sub_index_reg&(MASK_BIT(XRDY)))); Set_codec_cs_low(); MCBSP1_DXR1=codec_buf[i]; While(!(*Port_sub_index_reg&(MASK_BIT(XRDY)))); Set_codec_cs_high();} int ifilt,band,ns,n; { int k; double omega,lamda,epslon,f1,fh; double d[5],c[5]; void chebyi(),chebyii(),bwtf(); double coshl(),warp(),bpsub(),omin(); void fblt(); if((band==1)||(band==4)) fl=f1; if((band==2)||(band==3)) fl=f2; if(band<=3) fh=f3; if(band==4) fh=f4; if(ifilt<3) { switch(band) { case1: case2: { omega=warp(f2)/warp(f1); break; } case3: { omega=omin(bpsub(warp(f1),fh,fl),bpsub(warp(f4),fh,fl)); break; } case4: {omega=omin(1.0/bpsub(warp(f2),fh,fl),1.0/bpsub(warp(f3),fh,fl));} } lamda=pow(10.0,(db/20.0)); epslon=lamda/cosh(2*ns*coshl(omega)); } for(k=0;k{ switch(ifilt) { case1: { chebyi(2*ns,k,4,epslon,d,c); break; } case2: { chebyii(2*ns,k,4,omega,lamda,d,c); break; } case3: { bwtf(2*ns,k,4,d,c); break; } } fblt(d,c,n,band,fl,fh,&b[k*(n+1)+0],&a[(n+1)+0]); } } static double coshl(x) double x; { double z; z=log(x+sqrt(x*x-1.0)); return(z); } static double warp(f) double f; { double pi,z; pi=4.0*atan(1.0); z=tan(pi*f); return(z); } static double bpsub(om,fh,fl) double om,fh,fl; { double z; z=(om*om-warp(fh)*warp(fl))/((warp(fh)-warp(fl))*om); return(z); } static double omin(om1,om2) double om1,om2; { double z,z1,z2; z1=fabs(om1); z2=fabs(om2); z=(z1return(z); } static void bwtf(ln,k,n,d,c) int ln,k,n; double d[],c[]; { int i; double pi,tmp; pi=4.0*atan(1.0); d[0]=1.0; c[0]=1.0; for(i=1;i<=n;i++) { d[i]=0.0; c[i]=0.0; } tmp=(k+1)-(ln+1.0)/2.0; if(tmp==0,0) {c[1]=1.0;} else { c[1]=-2.0*cos((2*(k+1)+ln-1)*pi/(2*ln)); c[2]=1.0; } } static void chebyi(ln,k,n,ep,d,c) double d[],c[],ep; int ln,k,n; {int i; double pi,gam,omega,sigma; pi=4.0*atan(1.0); gam=pow(((1.0+sqrt(1.0+ep*ep))/ep),1.0/ln); sigma=0.5*(1.0/gam-gam)*sin((2*(k+1)-1)*pi/(2*ln)); omega=0.5*(1.0/gam+gam)*cos((2*(k+1)-1)*pi/(2*ln)); for(i=0;i<=n;i++) { d[i]=0.0; c[i]=0.0; } if(((ln%2)==1)&&((k+1)==(ln+1)/2)) { d[0]=-sigma; c[0]=d[0]; c[1]=1.0; } else { c[0]=sigma*sigma+omega*omega; c[1]=-2.0*sigma; c[2]=1.0; d[0]=c[0]; if(((ln%2)==0)&&(k==0)) d[0]=d[0]/sqrt(1.0+ep*ep); } } static void chebyii(ln,k,n,ws,att,d,c) double d[],c[],ws,att; int ln,k,n; { int i; double pi,gam,alpha,beta,sigma,omega,scln,scld; pi=4.0*atan(1.0); gam=pow((att+sqrt(att*att-1.0)),1.0/ln); alpha=0.5*(1.0/gam-gam)*sin((2*(k+1)-1)*pi/(2*ln)); beta=0.5*(1.0/gam+gam)*cos((2*(k+1)-1)*pi/(2*ln)); sigma=ws*alpha/(alpha*alpha+beta*beta); omega=-1.0*ws*beta/(alpha*alpha+beta*beta); for(i=0;i<=n;i++) { d[i]=0.0; c[i]=0.0; } if(((ln%2)==1)&&((k+1)==(ln+1)/2)) { d[0]=-1.0*sigma; c[0]=d[0]; c[1]=1.0; } else { scln=sigma*sigma+omega*omega; scld=pow((ws/cos((2*(k+1)-1)*pi/(2*ln))),2); d[0]=scln*scld; d[2]=scln; c[0]=d[0]; c[1]=-2.0*sigma*scld; c[2]=scld; } } #inlcude static void fblt(d,c,n,band,fln,fhn,b,a) int n,band; double fln,fhn,d[],d[],b[],a[]; { int i,k,m,n1,n2,ls; double pi,w,w0,w1,w2,tmp ,tmpd,tmpc,*work; double combin(); void bilinear(); pi=4.0*atan(1.0); w1=tan(pi*fln); for(i=n;i>=0;i--) { if((c[i]!=0.0)||(d[i]!=0.0)) break; } m=i; switch(band) { case1: case2: { n2=m; n1=n2+1; if(band==2) { for(i=0;i<=m/2;i++) { tmp=d[i]; d[i]=d[m-i]; d[m-i]=tmp; tmp=c[i]; c[i]c[m-i]; c[m-i]=tmp; } } for(i=0;i<=m;i++) { d[i]=d[i]/pow(w1,i); c[i]=c[i]/pow(w1,i); } break; } case3: case4: { n2=2*m; n1=n2+1; work=malloc(n1*n1*sizeof(double)); w2=tan(pi*fhn); w=w2-w1; w0=w1*w2; if(band==4) { for(i=0;i<=m/2;i++) { tmp=d[i]; d[i]=d[m-i]; d[m-i]=tmp; tmp=c[i]; c[i]=c[m-i]; c[m-i]=tmp; } } for(i=0;i<=n2;i++) { work[0*n1+i]=0.0; work[1*n1+i]=0.0; } for(i=0;i<=m;i++) { tmpd=d[i]*pow(w,(m-i)); tmpc=c[i]*pow(w,(m-i)); for(k=0;k<=i;k++) { ls=m+i-2*k; tmp=combin(i,i)/(combin(k,k)*combin(i-k,i-k)); work[0*n1+ls]+=tmpd*pow(w0,k)*tmp; work[1*n1+ls]+=tmpc*pow(w0,k)*tmp; } } for(i=0;i<=n2;i++) { d[i]=work[0*n1+i]; c[i]=work[1*n1+i]; } free(work); } } bilinear(d,c,b,a,n); } static double combin(i1,i2) int i1,i2; { int i; double s; s=1.0; if(i2==0) return(s); for(i=i1;i>(i1-i2);i--) {s*=i;} return(s); } static void bilinear(d,c,b,a,n) int n; double d[],c[],b[],a[]; { int i,j,n1; double sum,atmp,scale,*temp; n1=n+1; temp=malloc(n1*n1*sizeof(double)); for(j=0;j<=n;j++) {temp[j*n1+0]=1.0;} sum=1.0; for(i=1;i<=n;i++) { sum=sum*(double)(n-i+1)/(double)i; temp[0*n1+i]=sum; } for(i=1;i<=n;i++) for(j=1;j<=n;j++) { temp[j*n1+i]=temp[(j-1)*n1+i]-temp[j*n1+i-1]-temp[(j-1)*n1+i-1]; } for(i=n;i>=0;i--) { b[i]=0.0; atmp=0.0; for(j=0;j<=n;j++) { b[i]=b[i]+temp[j*n1+i]*d[j]; atmp=atmp+temp[j*n1+i]*c[j]; } scale=atmp; if(i!=0) a[i]=atmp; } for(i=0;i<=n;i++) { b[i]=b[i]/scale; a[i]=a[i]/scale; } a[0]=1.0; free(temp); }