In case Flash no longer exists; a copy of this site is included in the Flashpoint archive's "ultimate" collection.

Dead Code Preservation :: Archived AS3 works from wonderfl.net

forked from: Real-time Sound Spectrogram Visualizer

Real-time Sound Visualizer

@author Kosuke Suzuki
/**
 * Copyright rok.cancer ( http://wonderfl.net/user/rok.cancer )
 * MIT License ( http://www.opensource.org/licenses/mit-license.php )
 * Downloaded from: http://wonderfl.net/c/dy9O
 */

// forked from potato-attack's Real-time Sound Spectrogram Visualizer
package
{
  import adobe.utils.CustomActions;
  import flash.display.Bitmap;
  import flash.display.BitmapData;
  import flash.display.BlendMode;
  import flash.display.Sprite;
  import flash.display.StageAlign;
  import flash.display.StageQuality;
  import flash.display.StageScaleMode;
  import flash.events.Event;
  import flash.events.SampleDataEvent;
  import flash.geom.Point;
  import flash.geom.Rectangle;
  import flash.media.Microphone;
  import flash.media.Sound;
  import flash.media.SoundChannel;
  import flash.utils.ByteArray;
  import flash.utils.getTimer;

  /**
   * Real-time Sound Visualizer
   *
   * @author Kosuke Suzuki
   */
  public class Main extends Sprite
    {
      //----------------------------------------
      //CLASS CONSTANTS

      private const ZEROS:Point = new Point();




      //----------------------------------------
      //VARIABLES

      /**
       * マイクロフォン
       */
      private var _mic:Microphone;

      /**
       * 録音用
       */
      private var _records:Vector.<Number>;

      /**
       * 再生用
       */
      private var _sound:Sound;
      private var _soundChannel:SoundChannel;

      /**
       * 録音、再生ヘッダ
       */
      private var _position:int;

      /**
       * 描画用
       */
      private var _canvas:BitmapData;
      private var _slit:BitmapData;
      private var _over:Sprite;


      private var _pitch:PitchShifter;
      private var _fftFrameSize:int = 4096;

      private var windowValues        :Vector.<Number>;

      private var _channels:int = 60; // the number of pitch
      private var _constantQfilterbank:Vector.<Vector.<Number>>;
      private var _fl:Number = 50.0;  // The lowest frequency
      private var _Q:Number  = 60.0;  // Q value
      private var _fc:Number;         // center frequency
      //----------------------------------------
      //STAGE INSTANCES





      //----------------------------------------
      //METHODS

      /**
       * コンストラクタ
       */
      public function Main():void
        {
          // prepare constant-Q filterbank
          var _constantQfilterbank:Vector.<Vector.<Number>> = new Vector.<Vector.<Number>>();
          for (var m:int = 0; m < _channels; ++m)
            {
              _constantQfilterbank.push(new Vector.<Number>(_fftFrameSize/2));
              for (var n:int = 0; n< _fftFrameSize / 2; ++n)
                {
                  _fc = _fl * Math.pow(2, (m-1)/12);
                  _constantQfilterbank[m][n] = Math.exp( - (n-_fc)*(n-_fc) * _Q*_Q / 2 / _fc / _fc );
                }
            }
          
          //Wonderfl.disable_capture();

          stage.align = StageAlign.TOP_LEFT;
          stage.scaleMode = StageScaleMode.NO_SCALE;
          stage.quality = StageQuality.LOW;

          _mic = Microphone.getMicrophone();
          _mic.rate = 44;
          _mic.setSilenceLevel(0);
          _mic.setUseEchoSuppression(true);

          _sound = new Sound();

          _canvas = new BitmapData(stage.stageWidth, stage.stageHeight, false, 0x0000ff);
          _slit = new BitmapData(1, stage.stageHeight, false, 0x000000);
          addChild( new Bitmap(_canvas) );

          _over = addChild( new Sprite() ) as Sprite;
          _over.graphics.beginFill(0x0);
          _over.graphics.drawRect(0, 0, 100, 100);
          _over.graphics.endFill();

          _over.blendMode = BlendMode.INVERT; // invert the background color


          windowValues            = new Vector.<Number>(_fftFrameSize);
          //windowValuesFactored    = new Vector.<Number>(fftFrameSize);

          //var invFftFrameSize2:Number = 2.0 / (fftFrameSize2 * osamp
          var PI:Number        = Math.PI;
          var invFftFrameSizePI2:Number = PI * 2 / _fftFrameSize;
          for (var k:int = 0, t:Number = 0.0; k < _fftFrameSize; ++k, t += invFftFrameSizePI2)
            {
              var window: Number = -.5 * Math.cos(t) + .5;
              windowValues[k] = window;
              //windowValuesFactored[k] = window * invFftFrameSize2;
            }
          _pitch = new PitchShifter(_fftFrameSize, 4, 44100); // frame shift = 2048 / 4

          _startRecord();
        }

      /**
       * 録音開始
       */
      private function _startRecord():void
        {
          trace("_startRecord");

          _over.width = 0;

          _position = 0;
          //_records = new Vector.<Number>(44100 * 5);
          _records = new Vector.<Number>(_fftFrameSize);
          _mic.addEventListener(SampleDataEvent.SAMPLE_DATA, _micSampleDataHandler);
        }

      private var gFFTworksp    :Vector.<Number>; // FFT spectrum stack

      private var _magn:Number, _real:Number, _imag:Number; // variables for FFT
      private var _spectrogram_power : Vector.<Number>;
      private var _spectrum_temppower : Vector.<Number>;

      private var _canvasPosition:int = 0; // Position on Canvas

      private function _micSampleDataHandler(e:SampleDataEvent):void
        {
          var sw:int = stage.stageWidth;
          var sh:int = stage.stageHeight;

          var datas:ByteArray = e.data;
          gFFTworksp      = new Vector.<Number>(2 * _fftFrameSize + 2, true);

          while (datas.bytesAvailable)
            {
              var data:Number = _records[_position] = datas.readFloat();
              
              if (++_position == _records.length)
                {
                  _canvasPosition = _canvasPosition + 1;
                  
                  for (var k:int = 0, n:int = 1; k < _fftFrameSize; ++k, ++n)
                    {
                      gFFTworksp[n] = _records[k] * windowValues[k];
                      gFFTworksp[++n] = 0.0;
                    }
                  
                  _pitch.realft(gFFTworksp, _fftFrameSize, -1);
                  _position = 0;
                  
                  //
                  //                  private var _spectrogram_power : Vector.<Number>;
                  _spectrum_temppower = new Vector.<Number> ( _fftFrameSize / 2 );
                  
                  for (k = 0; k <= _fftFrameSize/2; ++k)
                    {
                      /* de-interlace FFT buffer */
                      _real = gFFTworksp[n = 1 + (k << 1)];
                      _imag = gFFTworksp[n + 1];
                      /* compute magnitude and phase */
                      _magn = 2.0 * Math.sqrt(_real * _real + _imag * _imag);

                      _spectrum_temppower[k] = _magn;
                    }

                  /*                    _mic.removeEventListener(SampleDataEvent.SAMPLE_DATA, _micSampleDataHandler);
                    _startSound();
                   */
                  //return;

                  // plot spectrogram
                  var x:Number = _canvasPosition % sw;
                  _canvas.copyPixels(_slit, _slit.rect, new Point(x + 1, 0));
                  for (var y:Number = 1; y <= _fftFrameSize / 2; ++y)
                    {
                        
                      _canvas.setPixel(x, sh - y,  Number( "0x" + _Number2RGB( Math.log(_spectrum_temppower[y]) / Math.LN10 ).toString(16) ) );
                    
                    }
                    trace(Math.log(_spectrum_temppower[200]) / Math.LN10);
                    
                }
              /*
              var x:Number = _position / _records.length * sw;
              _canvas.copyPixels(_slit, _slit.rect, new Point(x + 1, 0)); // erase previous wave at the next time
              
              var y:Number = sh * 0.5 - data * 1000;
              _canvas.setPixel(x, y, 0x00CC0000); // plot the point
               */
            }
        }

      /**
       *   Spectrum to Color RGB
       */
      private function _Number2RGB(_value:Number):uint
        {
          // var temp_value:Number;
          var max_value:Number =  3.0;
          var min_value:Number = -1.0;

          if (_value > max_value)
            {
              _value = max_value;
            }
          else if (_value < min_value)
            {
              _value = min_value;
            }

          // [-2,2] -> [0,1]
          _value = (_value - min_value) / (max_value - min_value);
          
          var h:Number = 240.0 * ( 1 - _value % 1 ) ; // normalize between 0 and 240
          
          var s:Number = 1.0;
          var v:Number = 1.0;
          
          var rgb:uint = 0;
          
          var hi:uint = Math.floor(h / 60.0) % 6;
          var f:Number = h / 60.0 - hi;
          var vv:uint = Math.round(255 * v);
          var pp:uint = Math.round(255 * v * ( 1 - s ));
          var qq:uint = Math.round(255 * v * ( 1 - f * s ));
          var tt:uint = Math.round(255 * v * ( 1 - (1 - f) * s ));
          if ( vv > 255 ) vv = 255;
          if ( pp > 255 ) pp = 255;
          if ( qq > 255 ) qq = 255;
          if ( tt > 255 ) tt = 255;
          switch (hi) {
          case 0: rgb = (vv << 16) | (tt << 8) | pp; break;
          case 1: rgb = (qq << 16) | (vv << 8) | pp; break;
          case 2: rgb = (pp << 16) | (vv << 8) | tt; break;
          case 3: rgb = (pp << 16) | (qq << 8) | vv; break;
          case 4: rgb = (tt << 16) | (pp << 8) | vv; break;
          case 5: rgb = (vv << 16) | (pp << 8) | qq; break;
          }
          
          // var _RGBst:uint = Number("0x" + temp_value.toString(16));
          //return _RGBst;
          return rgb;
        }
        
      /**
       * 再生開始
       */
      private function _startSound():void
        {
          trace("_startSound");

          _position = 0;
          _sound.addEventListener(SampleDataEvent.SAMPLE_DATA, _soundSampleDataHandler);
          _soundChannel = _sound.play();
        }

      private function _soundSampleDataHandler(e:SampleDataEvent):void
        {
          _over.width = stage.stageWidth * (_position / _records.length);
          _over.height = stage.stageHeight;

          for (var i:int = 0; i < 2048; ++i)
            {
              var data:Number = _records[_position];

              e.data.writeFloat(data);
              e.data.writeFloat(data);

              if (++_position == _records.length)
                {
                  _sound.removeEventListener(SampleDataEvent.SAMPLE_DATA, _soundSampleDataHandler);
                  _startRecord();
                  return;
                }
            }
        }
    }
}

/****************************************************************************
 *
 * NAME: PitchShifter.as
 * VERSION: 1.0
 * HOME URL: http://iq12.com/
 * KNOWN BUGS: none
 *
 * SYNOPSIS: Routine for doing pitch shifting while maintaining
 * duration using the Short Time Fourier Transform.
 *
 * DESCRIPTION: The routine takes a pitchShift factor value which is between 0.5
 * (one octave down) and 2. (one octave up). A value of exactly 1 does not change
 * the pitch. numSampsToProcess tells the routine how many samples in indata[0...
 * numSampsToProcess-1] should be pitch shifted and moved to outdata[0 ...
 * numSampsToProcess-1]. The two buffers can be identical (ie. it can process the
 * data in-place). fftFrameSize defines the FFT frame size used for the
 * processing. Typical values are 1024, 2048 and 4096. It may be any value <=
 * MAX_FRAME_LENGTH but it MUST be a power of 2. osamp is the STFT
 * oversampling factor which also determines the overlap between adjacent STFT
 * frames. It should at least be 4 for moderate scaling ratios. A value of 32 is
 * recommended for best quality. sampleRate takes the sample rate for the signal
 * in unit Hz, ie. 44100 for 44.1 kHz audio. The data passed to the routine in
 * indata[] should be in the range [-1.0, 1.0), which is also the output range
 * for the data, make sure you scale the data accordingly (for 16bit signed integers
 * you would have to divide (and multiply) by 32768).
 *
 * COPYRIGHT 1999-2006 Stephan M. Bernsee <smb [AT] dspdimension [DOT] com>
 *
 *                         The Wide Open License (WOL)
 *
 * Permission to use, copy, modify, distribute and sell this software and its
 * documentation for any purpose is hereby granted without fee, provided that
 * the above copyright notice and this license appear in all source copies.
 * THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY OF
 * ANY KIND. See http://www.dspguru.com/wol.htm for more information.
 *
 *****************************************************************************/

/****************************************************************************
 *
 * This code was converted to AS3/FP10 by Arnaud Gatouillat <fu [AT] iq12 [DOT] com>
 * from C# code by Michael Knight ( madmik3 at gmail dot com. )
 * http://sites.google.com/site/mikescoderama/
 *
 *****************************************************************************/

/****************************************************************************
 *
 * The functions `realft' and `four1' are based on those in Press, W.H., et al.,
 * Numerical Recipes in C: the Art of Scientific Computing (Cambridge Univ. Press,
 * 1989;  2nd ed., 1992).
 *
 *****************************************************************************/

class PitchShifter
{
  private var gInFIFO        :Vector.<Number>;
  private var gOutFIFO    :Vector.<Number>;
  private var gFFTworksp    :Vector.<Number>;
  private var gLastPhase    :Vector.<Number>;
  private var gSumPhase    :Vector.<Number>;
  private var gOutputAccum:Vector.<Number>;
  private var gAnaFreq    :Vector.<Number>;
  private var gAnaMagn    :Vector.<Number>;
  private var gSynFreq    :Vector.<Number>;
  private var gSynMagn    :Vector.<Number>;

  private var freqPerBin:Number, expct:Number;
  private var gRover:int, inFifoLatency:int, stepSize:int, fftFrameSize2:int;

  private var fftFrameSize:int, osamp:int, sampleRate:Number;

  /* pre-computed values for speed */
  private var windowValues        :Vector.<Number>;
  private var windowValuesFactored:Vector.<Number>;
  private var invPI:Number, invFftFrameSizePI2:Number, osampPI2:Number, invOsampPI2FreqBin:Number;

  private var PI:Number        = Math.PI
    private var TWOPI:Number    = 2 * Math.PI

      public function PitchShifter(fftFrameSize:int, osamp:int, sampleRate:Number)
        {
          this.fftFrameSize    = fftFrameSize;
          this.osamp            = osamp;
          this.sampleRate        = sampleRate;

          gInFIFO            = new Vector.<Number>(fftFrameSize);
          gOutFIFO        = new Vector.<Number>(fftFrameSize, true);
          gFFTworksp        = new Vector.<Number>(2 * fftFrameSize + 2, true);
          gLastPhase        = new Vector.<Number>(fftFrameSize / 2 + 1, true);
          gSumPhase        = new Vector.<Number>(fftFrameSize / 2 + 1, true);
          gOutputAccum    = new Vector.<Number>(2 * fftFrameSize, true);
          gAnaFreq        = new Vector.<Number>(fftFrameSize, true);
          gAnaMagn        = new Vector.<Number>(fftFrameSize, true);
          gSynFreq        = new Vector.<Number>(fftFrameSize, true);
          gSynMagn        = new Vector.<Number>(fftFrameSize, true);

          /* set up some handy variables */
          fftFrameSize2= fftFrameSize / 2;
          stepSize = fftFrameSize / osamp;
          freqPerBin = sampleRate / Number(fftFrameSize);
          expct = 2.0 * PI * Number(stepSize) / Number(fftFrameSize);
          inFifoLatency = fftFrameSize - stepSize;

          invPI = 1 / PI;
          invFftFrameSizePI2 = PI * 2 / fftFrameSize;
          osampPI2 = osamp / ( 2 * PI );
          invOsampPI2FreqBin = 1 / ( freqPerBin * osampPI2);

          windowValues            = new Vector.<Number>(fftFrameSize);
          windowValuesFactored    = new Vector.<Number>(fftFrameSize);

          var invFftFrameSize2:Number = 2.0 / (fftFrameSize2 * osamp);
          for (var k:int = 0, t:Number = 0.0; k < fftFrameSize; ++k, t += invFftFrameSizePI2)
            {
              var window: Number = -.5 * Math.cos(t) + .5;
              windowValues[k] = window;
              windowValuesFactored[k] = window * invFftFrameSize2;
            }
        }

  public function pitchShift(pitchShift:Number, numSampsToProcess:int, indata:Vector.<Number>):void
    {
      var magn:Number, phase:Number, tmp:Number, window:Number, real:Number, imag:Number, t:Number;
      var i:int, k:int, qpd:int, index:int, n:int;

      var outdata:Vector.<Number> = indata;
      if (gRover == 0) gRover = inFifoLatency;

      /* main processing loop */
      for (i = 0; i < numSampsToProcess; ++i)
        {
          /* As long as we have not yet collected enough data just read in */
          gInFIFO[gRover] = indata[i];
          outdata[i] = gOutFIFO[gRover - inFifoLatency];
          ++gRover;

          /* now we have enough data for processing */
          if (gRover >= fftFrameSize)
            {
              gRover = inFifoLatency;

              /* do windowing and re,im interleave */
              for (k = 0, n = 1; k < fftFrameSize; ++k, ++n)
                {
                  gFFTworksp[n] = gInFIFO[k] * windowValues[k];
                  gFFTworksp[++n] = 0.0;
                }
              /* ***************** ANALYSIS ******************* */
              /* do transform */
              realft(gFFTworksp, fftFrameSize, -1);
              /* this is the analysis step */
              for (k = 0; k <= fftFrameSize2; ++k)
                {
                  /* de-interlace FFT buffer */
                  real = gFFTworksp[n = 1 + (k << 1)];
                  imag = gFFTworksp[n + 1];

                  /* compute magnitude and phase */
                  magn = 2.0 * Math.sqrt(real * real + imag * imag);
                  phase = Math.atan2(imag, real);

                  /* compute phase difference */
                  tmp = phase - gLastPhase[k];
                  gLastPhase[k] = phase;

                  /* subtract expected phase difference */
                  tmp -= k * expct;

                  /* map delta phase into +/- Pi interval */
                  qpd = int(tmp * invPI);
                  if (qpd >= 0)    qpd += qpd & 1;
                  else            qpd -= qpd & 1;
                  tmp -= PI * Number(qpd);

                  /* get deviation from bin frequency from the +/- Pi interval */
                  tmp *= osampPI2;

                  /* compute the k-th partials' true frequency */
                  tmp = (k + tmp) * freqPerBin;

                  /* store magnitude and true frequency in analysis arrays */
                  gAnaMagn[k] = magn;
                  gAnaFreq[k] = tmp;

                }
              /* ***************** PROCESSING ******************* */
              /* this does the actual pitch shifting */
              for (var zero:int = 0; zero < fftFrameSize; ++zero)
                {
                  gSynMagn[zero] = 0.0;
                  gSynFreq[zero] = 0.0;
                }

              for (k = 0, n = pitchShift > 1.0 ? int(fftFrameSize2 / pitchShift) : fftFrameSize2; k <= n; ++k)
                {
                  index = int(k * pitchShift);
                  gSynMagn[index] += gAnaMagn[k];
                  gSynFreq[index] = gAnaFreq[k] * pitchShift;
                }
              /* ***************** SYNTHESIS ******************* */
              /* this is the synthesis step */
              for (k = 0; k <= fftFrameSize2; ++k)
                {
                  /* get magnitude and true frequency from synthesis arrays */
                  magn = gSynMagn[k];

                  /* subtract bin mid frequency */
                  /* get bin deviation from freq deviation */
                  /* take osamp into account */
                  /* add the overlap phase advance back in */
                  /* accumulate delta phase to get bin phase */
                  phase = (gSumPhase[k] += (gSynFreq[k] - Number(k) * freqPerBin) * invOsampPI2FreqBin + Number(k) * expct);

                  /* get real and imag part and re-interleave */
                  gFFTworksp[n = 1 + (k << 1)] = magn * Math.cos(phase);
                  gFFTworksp[n + 1] = magn * Math.sin(phase);
                }

              /* zero negative frequencies */
              for (k = fftFrameSize + 3, n = 1 + (fftFrameSize << 1); k < n; ++k)
                {
                  gFFTworksp[k] = 0.0;
                }
              /* do inverse transform */
              realft(gFFTworksp, fftFrameSize, 1);

              /* do windowing and add to output accumulator */
              for (k = 0, n = 1; k < fftFrameSize; ++k, ++n, ++n)
                {
                  gOutputAccum[k] += windowValuesFactored[k] * gFFTworksp[n];
                }
              for (k = 0; k < stepSize; ++k)
                {
                  gOutFIFO[k] = gOutputAccum[k];
                }

              //memmove(gOutputAccum, gOutputAccum + stepSize, fftFrameSize * sizeof(Number));
              /* shift accumulator */
              /* move input FIFO */
              for (k = 0, n = stepSize; k < inFifoLatency; ++k, ++n)
                {
                  gOutputAccum[k] = gOutputAccum[n];
                  gInFIFO[k] = gInFIFO[n];
                }
              for ( ;  k < fftFrameSize; ++k, ++n)
                {
                  gOutputAccum[k] = gOutputAccum[n];
                }
            }
        }
    }

  //private function realft( data:Vector.<Number>, n:int, isign:int ):void
  public function realft( data:Vector.<Number>, n:int, isign:int ):void
    {
      var i:int, i1:int, i2:int, i3:int, i4:int, n2p3:int;
      var c1:Number = 0.5, c2:Number, h1r:Number, h1i:Number, h2r:Number, h2i:Number;
      var wr:Number, wi:Number, wpr:Number, wpi:Number, wtemp:Number, theta:Number;

      theta = PI/n;
      if (isign == 1)
        {
          c2 = -0.5;
          four1(data, n, 1);
        }
      else
        {
          c2 = 0.5;
          theta = -theta;
        }
      wtemp = Math.sin(0.5 * theta);
      wpr = -2.0 * wtemp * wtemp;
      wpi = Math.sin(theta);
      wr = 1.0 + wpr;
      wi = wpi;
      n2p3 = 2 * n + 3;
      for (i = 2; i <= n / 2; ++i)
        {
          i4 = 1 + (i3 = n2p3 - (i2 = 1 + ( i1 = i + i - 1)));
          h1r =  c1 * (data[i1] + data[i3]);
          h1i =  c1 * (data[i2] - data[i4]);
          h2r = -c2 * (data[i2] + data[i4]);
          h2i =  c2 * (data[i1] - data[i3]);
          data[i1] =    h1r + wr * h2r - wi * h2i;
          data[i2] =    h1i + wr * h2i + wi * h2r;
          data[i3] =    h1r - wr * h2r + wi * h2i;
          data[i4] = -h1i + wr * h2i + wi * h2r;
          wr = (wtemp = wr) * wpr - wi * wpi + wr;
          wi = wi * wpr + wtemp * wpi + wi;
        }
      if (isign == 1)
        {
          data[1] = (h1r = data[1]) + data[2];
          data[2] = h1r - data[2];
        }
      else
        {
          data[1] = c1 * ((h1r = data[1]) + data[2]);
          data[2] = c1 * (h1r - data[2]);
          four1(data, n, -1);
          data=data;
        }
    }

  private function four1(data:Vector.<Number>, nn:int, isign:int):void
    {
      var n:int, mmax:int, m:int, j:int, istep:int, i:int;
      var wtemp:Number, wr:Number, wpr:Number, wpi:Number, wi:Number, theta:Number;
      var tempr:Number, tempi:Number;
      var j1:int, i1:int;
      n = nn << 1;
      j = 1;
      for (i = 1; i < n; i += 2)
        {
          if (j > i)
            {
              j1 = j + 1;
              i1 = i + 1;
              tempr = data[j];    data[j] = data[i];        data[i] = tempr;
              tempr = data[j1];    data[j1] = data[i1];    data[i1] = tempr;
            }
          m = n >> 1;
          while (m >= 2 && j > m)
            {
              j -= m;
              m >>= 1;
            }
          j += m;
        }
      mmax = 2;
      while (n > mmax)
        {
          istep = 2 * mmax;
          theta = TWOPI / (isign * mmax);
          wtemp = Math.sin(0.5 * theta);
          wpr = -2.0 * wtemp * wtemp;
          wpi = Math.sin(theta);
          wr = 1.0;
          wi = 0.0;
          for (m = 1; m < mmax; m += 2)
            {
              for (i = m; i <= n; i += istep)
                {
                  i1 = i +1;
                  j1 = 1+ (j = i + mmax);
                  tempr = wr*data[j]     - wi*data[j1];
                  tempi = wr*data[j1]  + wi*data[j];
                  data[j]     = data[i]     - tempr;
                  data[j1]    = data[i1]    - tempi;
                  data[i]  += tempr;
                  data[i1] += tempi;
                }
              wr = (wtemp = wr) * wpr - wi * wpi + wr;
              wi = wi * wpr + wtemp * wpi + wi;
            }
          mmax = istep;
        }
    }

}