Forums

Produce Audio Sound via Inverse FFT on Spectral Data

Started by louis June 13, 2005
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
louis wrote:
> I am > sure the numerical recipes code is reliable
I 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? )
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
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
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
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
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
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
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
>... i thought nobody was coding their own FFT no more. Do like >anybody, use FFTW, you'll never code anything as good as that
good 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