用C语言实现谱减法

该例程非常详细的实现了语音信号处理中的谱减法，具有很强的实践意义，希望可以帮到大家

/*---------------------------------------------------------------------
 * @（#）cfftall.c Version 1.1 by IKEDA Mikio 91/03/11
 * Comblex fft basic function
 * --------------------
 * cfftall（m0 x ainv）
 * --------------------
 * arguments
 *   int   m0   : 2**m0 = length
 *   double x[]  : complex signal --> spectrum
 *   double ainv : inverse switch （inverse = -1.0）
 *---------------------------------------------------------------------
 */

#include 
#include “ana.h“

void cfftall（int m0 double *x double ainv）
{
    int    i j lm li k lmx lmx2 np lix;
    double  temp1 temp2;
    double  c s csave sstep cstep;
    double  c0 s0 c1 s1;

    lmx = 1 << m0;

    csave = PI * 2.0 / （double）lmx;
    cstep = cos（csave）;
    sstep = sin（csave）;

    lmx += lmx;
    np   = lmx;

/*----- fft butterfly numeration */
    for （i = 3; i <= m0; ++i） {
	lix = lmx;
	lmx >>= 1;
	lmx2 = lmx >> 1;
	c = cstep;
	s = sstep;
	s0 = ainv * s;
	c1 = -s;
	s1 = ainv * c;
	for （li = 0; li < np; li += lix ） {
	    j = li;
	    k = j + lmx;
	    temp1  = x[j] - x[k];
	    x[j]  += x[k];
	    x[k]   = temp1;
	    temp2  = x[++j] - x[++k];
	    x[j]  += x[k];
	    x[k]   = temp2;

	    temp1  = x[++j] - x[++k];
	    x[j]  += x[k];
	    temp2  = x[++j] - x[++k];
	    x[j]  += x[k];
	    x[k-1] = c * temp1 - s0 * temp2;
	    x[k]   = s0 * temp1 + c * temp2;

	    j = li + lmx2;
	    k = j + lmx;
	    temp1  = x[j] - x[k];
	    x[j]  += x[k];
	    temp2  = x[++j] - x[++k];
	    x[j]  += x[k];
	    x[k-1] = -ainv * temp2;
	    x[k]   =  ainv * temp1;

	    temp1  = x[++j] - x[++k];
	    x[j]  += x[k];
	    temp2  = x[++j] - x[++k];
	    x[j]  += x[k];
	    x[k-1] = c1 * temp1 - s1 * temp2;
	    x[k]   = s1 * temp1 + c1 * temp2;

	}
	for （lm = 4; lm < lmx2; lm += 2） {
	    csave = c;
	    c = cstep * c - sstep * s;
            s = sstep * csave + cstep * s;

	    s0 = ainv * s;
	    c1 = -s;
	    s1 = ainv * c;

	    for （li = 0; li < np; li += lix ） {
		j = li + lm;
		k = j + lmx;
		temp1  = x[j] - x[k];
	        x[j]  += x[k];
		temp2  = x[++j] - x[++k];
		x[j]  += x[k];
	        x[k-1] = c * temp1 - s0 * temp2;
	        x[k]   = s0 * temp1 + c * temp2;

		j = li + lm + lmx2;
		k = j + lmx;
		temp1  = x[j] - x[k];
	        x[j]  += x[k];
		temp2  = x[++j] - x[++k];
		x[j]  += x[k];
	        x[k-1] = c1 * temp1 - s1 * temp2;
	        x[k]   = s1 * temp1 + c1 * temp2;
	    }
	}
	csave = cstep;
	cstep = 2.0 * cstep * cstep - 1.0;
	sstep = 2.0 * sstep * csave;
    }
    if （m0 >= 2）
	for （li = 0; li < np; li += 8） {
	    j = li;
	    k = j + 4;
	    temp1 = x[j] - x[k];
	    x[j] += x[k];
	    temp2 = x[++j] - x[++k];
	    x[j] += x[k];
	    x[k-1] = temp1;
	    x[k]   = temp2;
	    temp1  = x[++j] - x[++k];
	    x[j]  += x[k];
	    temp2  = x[++j] - x[++k];
	    x[j]  += x[k];
	    x[k-1] = -ainv * temp2;
	    x[k]   =  ainv * temp1;
	}
    for （li = 0; li < np; li += 4） {
	j = li;
	k = j + 2;
	temp1 = x[j] - x[k];
	x[j]  += x[k];
	x[k]   = temp1;
	temp2  

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2013-10-27 14:33  ss\
     文件       85762  2009-10-16 09:43  ss\5.pcm
     文件       85806  2009-10-13 12:22  ss\5.wav
     文件       85568  2009-10-16 09:49  ss\5out-1.pcm
     文件       85568  2009-10-16 09:36  ss\5out.pcm
     文件         398  1997-02-06 15:05  ss\Makefile
     文件          52  2000-06-01 20:20  ss\README
     文件        1226  1997-02-06 14:02  ss\ana.h
     文件        3521  2013-10-27 14:33  ss\cfftall.c
     文件         898  1997-02-06 13:56  ss\getfirst.c
     文件        1282  1997-02-06 13:56  ss\irfft.c
     文件         632  1997-02-06 13:56  ss\multirr.c
     文件      142784  2009-10-16 11:10  ss\out.pcm
     文件         945  1997-02-06 13:56  ss\rdframe.c
     文件        1680  1997-02-06 13:57  ss\rfft.c
     文件        6242  2009-10-21 15:56  ss\spsub.c
     文件      237760  2009-10-16 12:21  ss\tst1-yanhua-out.pcm
     文件      237804  2009-10-16 11:48  ss\tst1-yanhua-out.wav
     文件      237760  2009-10-16 11:36  ss\tst1-yanhua-out1.pcm
     文件       91776  2009-10-16 11:05  ss\tst1out-1.pcm
     文件         635  1997-02-06 14:00  ss\xd_realloc.c