void PulseWidth::processVariable(int numFrames, int offsetFrame, long framesProcessed) {
    float *output = outputBuffer(0)->offset(offsetFrame);
    float *table = mWaveTable.offset(0);
    float *delayBuffer = mDelayBuffer.offset(0);

    
    float sr = sampleRate();    

    float nyq = 0.5; // number of harmonics in the band limited pulse
    int n;
    float nyqHz = sr * nyq;
        
    float sqr;

    float *freq = mFrequency->buffer()->offset(offsetFrame);
    float fqc = fabs( freq[0] );
    if( fqc <= 0.001 ) fqc = 0.001;
    if( fqc >= nyqHz ) fqc = nyqHz;
        
    float pw = mPulseWidth->value();

    // harmonics
    int knh = (int)(sr * nyq / fqc);
    if((n = (int)knh) <= 0) n = 1;    
    long tnp1 = n + n + 1;   // calc 2n + 1
    float over2n = 0.5 / n;


    float pulse = 0.0;

    float out1;
    
    int downPhase;
    
    float scale = over2n, denom, num;

    // for interpolation
    long  v1, v2;    
    float fv1;
     
    long inc = (long)(fqc * mSicvt2);
    
//    printf("framesProcessed=%d\n", framesProcessed);
    
    for(int i = 0; i < numFrames; i++) {
        
//        printf("%f ", fqc);
        
        downPhase = mPhase >> mWaveTable.lobits;
        denom = *(table + downPhase);

        if(denom > 0.00001 || denom < -0.00001) {
            num = *(table + (downPhase * tnp1 & mWaveTable.lenmask));
            pulse = (num / denom - 1.0) * scale;
        } else {
            pulse = 1.0;
        }

        mPhase += inc;
        mPhase &= PHMASK;

        fqc = fabs( freq[i] );
        if( fqc <= 0.001 ) fqc = 0.001;
        if( fqc >= nyqHz ) fqc = nyqHz;
            
        // harmonics
        knh = (int)(sr * nyq / fqc);
        if((n = (int)knh) <= 0) n = 1;    
        tnp1 = n + n + 1;   // calc 2n + 1
        over2n = 0.5 / n;
            
                    
        inc  = (long)(fqc* mSicvt2);

        // VDelay inserted here
        delayBuffer[mDelayIndex] = pulse;
        fv1 = (float) mDelayIndex - sr * pw / fqc;

        v1 = (long) fv1;
        if(fv1 < 0.0) v1--;

        fv1 -= (float) v1;

        // Make sure Inside the buffer
        while (v1 >= mMaxDelaySamples) v1 -= mMaxDelaySamples;        
        while (v1 < 0L) v1 += mMaxDelaySamples;

        // Find next sample for interpolation
        v2 = (v1 < (mMaxDelaySamples - 1L) ? v1 + 1L : 0L);
        out1 = delayBuffer[v1] + fv1 * (delayBuffer[v2] - delayBuffer[v1]);

        if (++mDelayIndex == mMaxDelaySamples) mDelayIndex = 0;
        // End of VDelay


        sqr  = pulse - out1 + mLeaky * mYnm1;
        mYnm1 = sqr;
        
        // the 1.9 normalizes (but does initially clip)
        float o = (sqr + pw - 0.5); // * 1.9;
        output[i]  = o;
    }
    
}