본문 바로가기
소프트웨어개발

[ C ] 소스 IPC_spectrum.c

by inkul 2014. 8. 23.
  1. /*Introducing Speech and Language Processing  
  2. http://www.islp.org.uk*/   
  3.    
  4.    
  5. //lpc_spectrum.c   
  6. ===============================================   
  7. #include <stdio.h>   
  8. #include <math.h>   
  9. #include <stdlib.h>   
  10. #include "slputils.c"   
  11. #include "four1.c"   
  12.    
  13. /* LPC_SPECTRUM.C                                                                */   
  14. /* Reads a signal from a disk file into a variable, x_in.                        */   
  15. /* Computes the LPC coefficients of sample_number, pads it with zeroes to 512        */   
  16. /* and writes the resulting signal to outfile.                                */   
  17.    
  18. int main(int argc, char *argv[])   
  19. {   
  20.         short int *x_in, *signal_in();   
  21.         char *infile, *sampleno, *endcp;   
  22.         int sample, i, *length=0;   
  23.         float frame[81], *xms, d[17], *x_in_f, logpsd[512], x_frame[512];   
  24.         void memcof();   
  25.         float fftdata[1024];                        /* 512 complex data points        */   
  26.         void four1();   
  27.    
  28.         if (argc != 3) {   
  29.                 printf("usage: lpc_spectrum input_file sample_number [> output_file]\n");   
  30.                 exit(1);   
  31.         }   
  32.    
  33.         infile = argv[1];   
  34.         sampleno = argv[2];   
  35.    
  36.         x_in = signal_in(infile,length);   
  37.         sample = (int) strtoul(sampleno,&endcp,10);   
  38.    
  39. /* Make a floating-point version of x_in, called x_in_f                        */   
  40.         x_in_f = (float *) calloc(*length,sizeof(float));   
  41.         if (!x_in_f) {   
  42.                 printf("Can't allocate space for x_in_f.\n");   
  43.                 exit(1);   
  44.         }   
  45.    
  46.         for (i=0;i<*length;i++) {   
  47.                 x_in_f[i] = ((float) x_in[i]);   
  48.         }   
  49.    
  50.         for (i=0;i<=79;i++) {   
  51.                 frame[i+1] = x_in_f[sample+i];   
  52.         }   
  53.    
  54.         memcof(frame,80,10,&xms,d);   
  55.    
  56.         for(i=0;i<251;i++) x_frame[i] = 0;    
  57.         for(i=0;i<10;i++) x_frame[i+248] = 15000.0*d[i+1];    
  58.         for(i=261;i<512;i++) x_frame[i] = 0;   
  59.    
  60.         for (i=0;i<=511;i++) {   
  61.                 fftdata[2*i] = (float) x_frame[i];   
  62.                 fftdata[2*i+1] = 0.0;   
  63.         }   
  64.    
  65.         four1(fftdata-1,512,1);   
  66.    
  67. /* In the log power spectral density, magnitudes are in dB                */   
  68. /* in steps of SampleRate/512 Hz (15.625Hz at 8000 samples/s).        */   
  69.    
  70.    printf("f (Hz)\tAmplitude (dB)\n");   
  71.    
  72.    for (i=0;i<=256;i++) {   
  73.       logpsd[i] = fftdata[2*i];   
  74.       logpsd[i] *= logpsd[i];   
  75.       logpsd[i] += SQR(fftdata[2*i+1]);   
  76.       logpsd[i] = -10*log10(logpsd[i]);   
  77.       printf("%d\t%.2f\n",(int) (i*15.625),logpsd[i]);   
  78.    }   
  79.    return 0;   
  80. }   
  81.    
  82.    
  83. four1.c   
  84. ===============================================   
  85. #include <math.h>   
  86. #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr   
  87.    
  88. /* four1.c                                                                */   
  89. /* Discrete Fourier Transform                                                */   
  90. /* This routine is from the book Numerical Recipes in C.                */   
  91. /* (Cambridge University Press), Copyright (C)  1987-1992 by            */   
  92. /* Numerical Recipes Software.  Used by permission.  Use of this        */   
  93. /* routine other than as an integral part of the software collection    */   
  94. /* provided in accompaniment with the book Introducing Speech and       */   
  95. /* Language Processing (Cambridge University Press) requires an         */   
  96. /* additional license from Numerical Recipes Software ([url]www.nr.com).[/url]     */   
  97. /* Further distribution in any form is prohibited.                      */   
  98.    
  99. void four1(float data[], unsigned long nn, int isign)   
  100. /* Replaces data[1..2*nn] by its discrete Fourier transform, if isign        */   
  101. /* is input as 1; or replaces data[1..2*nn] by nn times its inverse        */   
  102. /* discrete Fourier transform, if isign is input as -1. data is a         */   
  103. /* complex array of length nn or, equivalently, a real array of length        */   
  104. /* 2**nn. nn MUST be an integer power of 2 (this is not checked for!).        */   
  105. {   
  106.         unsigned long n,mmax,m,j,istep,i;   
  107.         double wtemp,wr,wpr,wpi,wi,theta;   
  108.         float tempr,tempi;   
  109.    
  110.         n=nn << 1;   
  111.         j=1;   
  112.         for(i=1;i<n;i+=2) {   
  113.                 if (j > i) {   
  114.                         SWAP(data[j],data[i]);   
  115.                         SWAP(data[j+1],data[i+1]);   
  116.                 }   
  117.                 m=n >> 1;   
  118.                 while (m >= 2 && j > m) {   
  119.                         j -= m;   
  120.                         m >>= 1;   
  121.                 }   
  122.                 j += m;   
  123.         }   
  124.            
  125.         mmax=2;   
  126.         while(n > mmax) {   
  127.                 istep=mmax << 1;   
  128.                 theta=isign*(6.28318503717959/mmax);   
  129.                 wtemp=sin(0.5*theta);   
  130.                 wpr = -2.0*wtemp*wtemp;   
  131.                 wpi=sin(theta);   
  132.                 wr=1.0;   
  133.                 wi=0.0;   
  134.                 for(m=1;m<mmax;m+=2) {   
  135.                         for (i=m;i<=n;i+=istep) {   
  136.                                 j=i+mmax;   
  137.                                 tempr=wr*data[j]-wi*data[j+1];   
  138.                                 tempi=wr*data[j+1]+wi*data[j];   
  139.                                 data[j]=data[i]-tempr;   
  140.                                 data[j+1]=data[i+1]-tempi;   
  141.                                 data[i] += tempr;   
  142.                                 data[i+1]+= tempi;   
  143.                         }   
  144.                         wr=(wtemp=wr)*wpr-wi*wpi+wr;   
  145.                         wi=wi*wpr+wtemp*wpi+wi;   
  146.                 }   
  147.                 mmax=istep;   
  148.         }   
  149. }   
  150.    
  151.    
  152.    
  153. slputils.c   
  154. ===============================================   
  155. #include <stdio.h>   
  156. #include <stdlib.h>   
  157. #include "nrutil.h"   
  158.    
  159. /* SLPUTILS.C                                                                        */   
  160.    
  161. short int *signal_in(char *infile, int *length)   
  162.    
  163. /* Reads a signal from a disk file into a variable, x_in.                        */   
  164. /* The length of the signal (in shorts) is returned in the variable *length.        */   
  165.    
  166. {   
  167.         FILE *fid;   
  168.         int j, status;   
  169.         short int input, *x_in;   
  170.    
  171. /* Open the file */   
  172.    
  173.         fid = fopen(infile,"rb");   
  174.         if(fid == NULL) {   
  175.                 printf("Error opening %s\n",infile);   
  176.                 exit(1);   
  177.         }   
  178.    
  179. /* First read in shorts, until the end of the file is reached, just to measure the file */   
  180.         j = 0;   
  181.         while(!feof(fid)) {   
  182.                 status = fread(&input,sizeof(short int),1,fid);   
  183.                 if(status != 1) {   
  184.                         break;   
  185.                 }   
  186.                 j++;   
  187.         }   
  188.         (*length) = j;   
  189.         fclose(fid);   
  190.    
  191. /* Allocate memory, and read in the file again        */   
  192.    
  193.         x_in = (short int *) calloc(*length,sizeof(short int));   
  194.    
  195. /* Open the file again */   
  196.    
  197.         fid = fopen(infile,"rb");   
  198.         if(fid == NULL) {   
  199.                 printf("Error re-opening %s\n",infile);   
  200.                 exit(1);   
  201.         }   
  202.    
  203.         for(j=0;j<*length;j++) {   
  204.                 status = fread(&input,sizeof(short int),1,fid);   
  205.                 x_in[j] = input;   
  206.                 if(status != 1) {   
  207.                         printf("Read error at sample %d\n",j);   
  208.                         exit(1);   
  209.                 }   
  210.         }   
  211.         *length = j;   
  212.         fclose(fid);   
  213.         return(x_in);   
  214. }   
  215.    
  216. void signal_out(int *length, short int *x_out, char *outfile)   
  217. {   
  218.         FILE *fid;   
  219.         int status;   
  220.    
  221.         fid = fopen(outfile,"wb");   
  222.         if(fid == NULL) {   
  223.                 printf("Can't open %s\n",outfile);   
  224.                 printf("It may be in use by another application.\n");   
  225.                 exit(1);   
  226.         }   
  227.         status = fwrite(x_out,sizeof(short int),*length,fid);   
  228.         if(status < *length) {   
  229.                 printf("Error writing %s\n",outfile);           
  230.                 exit(1);   
  231.         }   
  232.         fclose(fid);   
  233. }   
  234.    
  235.    
  236.    
  237. nrutil.h   
  238. ===============================================   
  239. /* nrutil.h                                                                */   
  240.    
  241. /* From W. H.Press, S. A. Teukolsky, W. T. Vettering and B. P. Flannery */   
  242. /* (1992) Numerical Recipes in C: The Art of Scientific Computing        */   
  243. /* (2nd edition). Cambridge University Press.                                */   
  244.    
  245. /* Non-copyright                                                         */   
  246.    
  247. /* pp. 941-2: This file only contains selected portions, however.        */   
  248.    
  249. static float sqrarg;   
  250. #define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)   
  251.    
  252. static float maxarg1,maxarg2;   
  253. #define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))   
  254.    
  255. void nrerror(char error_text[]);   
  256. float *vector(long nl, long nh);   
  257. void free_vector(float *v, long nl, long nh); 

댓글0