Hello, I am using C++ on a linux little endian machine, and am interfacing to OpenSound system to produce my audio sound. So my goal is to be able to reprduce the correct raw sound data corresponding to spectral data I already have (for example, I know the frequency versus power information for at least 128 bins, and over about 30 time points per second) I am using the four1 function found on the Numerical Recipes in C pages online (their FFT code is posted online). So before I do this, I am first trying to output a single frequency tone to make sure the IFFT is working - but low and behold it doesnt work. I am sure the numerical recipes code is reliable as many people have experience with it, so maybe I am setting it up wrong? Maybe I don't have the imaginary and real components set right? or maybe its even a type casting thing? Who knows! Below is the main code: (The last few lines outputs the raw sound data, and I later play this back using the correct channel and sampling rate info, however I don't hear a single frequency tone like I should) static unsigned char buf[30*1024]; static unsigned char large_buffer_audio[30*1024]; static unsigned char buf_audio[30*1024]; void four1(float [], unsigned long, int ); //the third parameter in four1 corresponds to forward or // backward fft, if set to -1 then do inverse fft int main() { const int fft_size = 1024; float f1[2048]; int i, j, k, t; FILE *fd_audio; for (t = 0; t<30; t++ ) { for (k=0; k<fft_size; k++) { f1[ 100 ] = 1000; f1[fft_size-100] = 1000.0; f1[k] = 0; } four1(f1-1, fft_size, -1); for (i=0; i< 2*fft_size; i++) { buf_audio[ i ] = (unsigned char)(f1[(2*i)]); } for ( i = (t*1024), j=0; j< 2*fft_size;i++,j++ ) { large_buffer_audio[ i ] = buf_audio[j]; large_buffer_audio[ i ] = buf_audio[j]; } } FILE *audiofile = fopen( "myfile_num_rec_fft", "wb" ); fwrite(large_buffer_audio,sizeof(unsigned char),sizeof(large_buffer_audio), audiofile ); fclose(fd_audio); return 0; } Thank you for your assistance!!! Its very appreciated! LD This message was sent using the Comp.DSP web interface on www.DSPRelated.com

# Produce Audio Sound via Inverse FFT on Spectral Data

Started by ●June 13, 2005

Reply by ●June 14, 20052005-06-14

louis wrote:> I am > sure the numerical recipes code is reliableI wouldn't be so quick to assume that. I'd recommend a little googling. I, for one, will never trust NR again. There are other choices. My kissfft library might be a good match for your needs. http://www.sourceforge.net/projects/kissfft It will probably be fast enough, but if you need faster and can tolerate GPL in your code -- FFTW is a great library (fftw.org). There are also C-calllable routines in GSL for FFT (also GPL). If I were doing your project, I might use Python. The Numeric module has FFT support built in. Maybe I don't have the> imaginary and real components set right? or maybe its even a type casting > thing? Who knows! Below is the main code: > (The last few lines outputs the raw sound data, and I later play this back > using the correct channel and sampling rate info, however I don't hear a > single frequency tone like I should) > > static unsigned char buf[30*1024]; > static unsigned char large_buffer_audio[30*1024]; > static unsigned char buf_audio[30*1024]; > > void four1(float [], unsigned long, int ); > //the third parameter in four1 corresponds to forward or > // backward fft, if set to -1 then do inverse fft > > int main() > { > const int fft_size = 1024; > float f1[2048]; > int i, j, k, t; > FILE *fd_audio; > > for (t = 0; t<30; t++ ) > { > for (k=0; k<fft_size; k++) > { > f1[ 100 ] = 1000; > f1[fft_size-100] = 1000.0; > f1[k] = 0;> } ^^^^^^^^^^ I don't think this is what you meant to do. Actually it looks like you never even initialize the latter half of your array.> > four1(f1-1, fft_size, -1);^^^^ why f1-1 ?? Oh I think I know -- is this b/c the NR people couldn't be bothered to actually convert their Fortran code from one-based indexing to zero-based? ( Are you starting to get the idea that I don't think much of Numerical Recipes? )

Reply by ●June 20, 20052005-06-20

Hi thanks very much for your response, Okay below I have corrected my previously posted code - Here I outputting a 1khz tone, using the Numerical Recipes four1 code. Note that I am using a 1024 point FFT, where I specify that my frequencies range from +22050 to -22050, so each frequency bin is 43 Hz. To output a 1 Khz tone, then I need to assign a positive amplitude to bin 1000/43 = 23 I noticed that the code below only works for when I have even numbered bins - so for example if I change the '23' to a '24' , then this will output a beautiful sounding 1khz tone. However if i leave it at bin 23, then there is some noise. I tried this for numerous bins, and the pattern holds - even bins work, odd ones don't. I feel like i'm really close - and I am probably not assigning something correctly? But I've looked at this a million times and don't know for the life of me what the problem could be! four1 really should work..!! Is a c++ version of an ifft that might work better?? Thank you in advance!! int main() { const int fft_size = 1024; float f1[2*fft_size]; float f2[2*fft_size]; float *ff1, *ff2; int i, j, k, t; ff1 = f1-1; ff2 = f2-1; for ( t=0; t<79; t++) { for( k=0; k<2*fft_size; k++) { f1[ k ] = 0.0; f2[ k ] = 0.0; f1[ (23*2) ]= 1000.0; //1000Hz/43Hz = bin 23 } four1( ff1, (unsigned long)1024, -1); four1( ff2, (unsigned long)1024, -1); for (i=0; i< fft_size; i++) { buf_audio[ 2*i ] = (short int)(f1[(2*i)]); // left channel buf_audio[ (2*i) + 1 ] = (short int)(f2[(2*i)]); // right channel } for ( i = (t*1024), j=0; j < fft_size ; i++, j++ ) { large_buffer_audio[ i ] = buf_audio[j]; large_buffer_audio[ i ] = buf_audio[j]; } } FILE *audiofile = fopen( "debug1", "wb" ); fwrite( large_buffer_audio, sizeof(short int), sizeof(large_buffer_audio), audiofile ); return 0; } This message was sent using the Comp.DSP web interface on www.DSPRelated.com

Reply by ●June 20, 20052005-06-20

OKAY so it turns out it actually works, it was a mistake in the way i was outputting the sound data anyways, now ive applied this to spectral data that i have, but it doesnt produce the right sound theoretically shouldnt this work? This message was sent using the Comp.DSP web interface on www.DSPRelated.com

Reply by ●June 21, 20052005-06-21

louis wrote:> theoretically shouldnt this work?>From your original post it seems that you only have the power(magnitude) information in spectral domain. In this case, no, theoretically it shouldn't work. Regards, Andor

Reply by ●June 21, 20052005-06-21

okay thank you again for your post, but can you explain why theoretically it can't work?? (i have taken a dsp course, but i feel my signal processing background is still weak) so just to clarify - when i assign a positive amplitude to a bin that corresponds to say frequency 1000hz, then do the ifft, i hear a clear sounding 1khz tone. Thne when i assign two frequency bins, say bins 3 and 4 (corresponding to frequencies 129 and 172 hz, since each bin is 43 hz) then i get some waveform. I used matlab to plot the summation of two sine waves 129 and 172 hz, and the waveform looks different. My plot doesnt look noisy or anything, but its not right when i compare it to matlab.. So i dont know what to make of all this Is it me? is it Numerical recipes?? My boss is insistent that it shoudl be possible :s (i am ready to have a nervous breakdown!!) okay thank you kindly again for your assistance!! This message was sent using the Comp.DSP web interface on www.DSPRelated.com

Reply by ●June 22, 20052005-06-22

louis wrote:> okay thank you again for your post, but can you explain why theoretically > it can't work??Because you need magnitude and phase, or equivalently real and imaginary part, to fully specify the DFT X of a time domain signal x. For example sin(t), -sin(t) and cos(t) all have the same magnitude but different phase. Regards, Andor

Reply by ●June 23, 20052005-06-23

louis allegedly wrote:> so just to clarify - when i assign a positive amplitude to a bin that > corresponds to say frequency 1000hz, then do the ifft, i hear a clear > sounding 1khz tone. Thne when i assign two frequency bins, say bins 3 and > 4 (corresponding to frequencies 129 and 172 hz, since each bin is 43 hz) > then i get some waveform. I used matlab to plot the summation of two sine > waves 129 and 172 hz, and the waveform looks different.Depending on whether you assign to the real or imaginary bins (in a complex IFFT, there are 2 bins for each frequency), you could get the sum of 2 sine waves, 2 cosine waves, or a sin+cos waveform (etc.), which will all look different. IMHO. YMMV. -- rhn Ron Nicholson

Reply by ●June 24, 20052005-06-24

damn, i thought nobody was coding their own FFT no more. Do like anybody, use FFTW, you'll never code anything as good as that louis wrote:> Hello, > > I am using C++ on a linux little endian machine, and am interfacing to > OpenSound system to produce my audio sound. So my goal is to be able to > reprduce the correct raw sound data corresponding to spectral data I > already have (for example, I know the frequency versus power information > for at least 128 bins, and over about 30 time points per second) > > I am using the four1 function found on the Numerical Recipes in C pages > online (their FFT code is posted online). > > So before I do this, I am first trying to output a single frequency tone > to make sure the IFFT is working - but low and behold it doesnt work. I am > sure the numerical recipes code is reliable as many people have experience > with it, so maybe I am setting it up wrong? Maybe I don't have the > imaginary and real components set right? or maybe its even a type casting > thing? Who knows! Below is the main code: > (The last few lines outputs the raw sound data, and I later play this back > using the correct channel and sampling rate info, however I don't hear a > single frequency tone like I should) > > static unsigned char buf[30*1024]; > static unsigned char large_buffer_audio[30*1024]; > static unsigned char buf_audio[30*1024]; > > void four1(float [], unsigned long, int ); > //the third parameter in four1 corresponds to forward or > // backward fft, if set to -1 then do inverse fft > > int main() > { > const int fft_size = 1024; > float f1[2048]; > int i, j, k, t; > FILE *fd_audio; > > for (t = 0; t<30; t++ ) > { > for (k=0; k<fft_size; k++) > { > f1[ 100 ] = 1000; > f1[fft_size-100] = 1000.0; > f1[k] = 0; > } > > four1(f1-1, fft_size, -1); > > for (i=0; i< 2*fft_size; i++) > { > buf_audio[ i ] = (unsigned char)(f1[(2*i)]); > } > > > for ( i = (t*1024), j=0; j< 2*fft_size;i++,j++ ) > { > large_buffer_audio[ i ] = buf_audio[j]; > large_buffer_audio[ i ] = buf_audio[j]; > } > } > > > FILE *audiofile = fopen( "myfile_num_rec_fft", "wb" ); > fwrite(large_buffer_audio,sizeof(unsigned > char),sizeof(large_buffer_audio), audiofile ); > > fclose(fd_audio); > return 0; > } > > > Thank you for your assistance!!! Its very appreciated! > LD > > > > > This message was sent using the Comp.DSP web interface on > www.DSPRelated.com

Reply by ●June 24, 20052005-06-24

>... i thought nobody was coding their own FFT no more. Do like >anybody, use FFTW, you'll never code anything as good as thatgood for what purpose? If you code your own FFT, you understand much better how it works (or doesn't work), and you own the copyright. If you design your own hardware (a custorm VLIW inside an FPGA for instance) you might be able to optimize for your specific constraints better than by just letting a compiler do its job on canned code. IMHO. YMMV. -- rhn Ron Nicholson