[ C ] 소스 IPC_spectrum.c

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);