Construction of Room Impulse Response - Short-Time-Fourier-Transform and Inverse-STFT - Matlab and NWaves (C#) Implementation

Question

Background: I am working on a Unity Project (C#) and implementing a Room Impulse Response synthesis. I do ray tracing to create energy histogram for different bands. With these histograms it is possible to synthesize a room impulse response.

My guide is following theoretical Reference: Physically based real-time auralization of interactive virtual environments (Schröder) 5.3.4 Construction of the Room Impulse Response (RT) https://publications.rwth-aachen.de/record/50580/files/3875.pdf

I also found a helpful Matlab Implementation which refers to the work above: https://de.mathworks.com/help/audio/ug/room-impulse-response-simulation-with-stochastic-ray-tracing.html#RoomImpulseWithRayTracingExample-8

I worked through the process step by step and so far it managed to implement this to my C# class. Now I am stuck at the point where I have to use the STFT and the ISTFT. For that reason I am using the NWaves (C#) package to work with (I)STFT: https://github.com/ar1st0crat/NWaves/wiki/Transforms#short-time-fourier-transform

To be more precise, I have (poisson distributed) dirac pulses in a 44100 Array (1sek IR 44100 samplerate). This time signal has to be transformed to frequency domain via STFT to filter it with a Raised Cosine Filter and transform it back to time domain.

In the matlab code its at the part Filter the Poisson sequence through the six bandpass filters. Without going into further details, the problem is that the Matlab Function sfft(signalVector) is returning a different result than NWaves. Right now I am using the Matlab Example Online Editor to evaluate its output (via Try this example on the website).

I don't understand how to (1. )interpret the results or (2.) replicate the calculation process of Matlab. It seems like Matlab is doing so much under the hood. With no experience in Matlab this is very hard to comprehend.

1.) I have isolated the problem to grasp it but Matlab and NWaves are returning different results despite the same input and I don't understand why. I am using the same STFT configuration for Matlab and NWaves.

Reduced Matlab Code

NFFT = 8192;
win = hann(882,"symmetric");
sfft = dsp.STFT(Window = win,OverlapLength=441,FFTLength=NFFT,FrequencyRange="onesided");
isfft = dsp.ISTFT(Window=win,OverlapLength=441,FrequencyRange="onesided");
x = linspace(1, 1, 441).'
X = sfft(x);

Matlab returns a 4097 Array with a real and imaginary part. I expected a length of 4097 due to the FFT Length (FFT Size / 2 + 1) As far I know the imaginary values are strictly mandatory to reconstruct the frequency domain signal back to time domain.

The reduced C# Code:

float[] testValues = new float[441];
for (int i = 0; i < 441; i++)
{
    testValues[i] = 1.0f;
}

int nFFT = 8192; // FFT Size defines the frequency resolution of the spectrogram. 
int hannWindowLength = 882; // The window length defines the time resolution of the spectrogram.
int hopSize = hannWindowLength / 2; // The Hopsize defines the overlap between two consecutive windows. The hopsize is usually chosen to be 50% of the window length. 
Stft stft = new Stft(hannWindowLength, hopSize, NWaves.Windows.WindowType.Hann, nFFT);

List < (float[], float[]) > stftFrameResult = stft.Direct(testValues);
float[] realSpectrum1 = stftFrameResult[0].Item1; // Length: 8192
float[] imagSpectrum1 = stftFrameResult[0].Item2; // Length: 8192

for (int j = 0; j < realSpectrum1.Length; j++)
{
    Debug.Log(realSpectrum1[j] + " " + imagSpectrum1[j]);
}

NWaves returns a complete different result. The spectrum has the length of the FFT (8192) and all values except the first one are different.

In NWaves I can not set the frequency range. Therefore I interprete its results that I have to use only the first half of the Array to 4097. But the results would be still different and I don't understand why.

When I change the FrequencyRange in Matlab to twosided I receive a result of length 8192. But the values still differs from NWaves.

2.) Even if the results would be the same. How is Matlab doing his sfft magic? In the linked example (Filter the Poisson sequence through the six bandpass filters)

 X = sfft(x);
 X = X.*RCF.';
 y((index-1)*frameLength+1:index*frameLength,:) = isfft(X);

Where RCF is a 6 x 4097 Matrix with only real values. X will become a 4097x6 Matrix with real and imaginary parts due to filtering with the RCF. Are these imaginary parts somehow included in the calculation of X.*RCF? Or are they only attached to the matrix but not touched?

I would appreciate any share of knowledge or hints. Thanks!

Answer 1

If someone will ever stumble on this, here is the creation of the dirac sequence + filtering it via the RCF using nWaves for the STFT

    private float[, ] CreateAndFilterDiracSequenceWithRCF(float roomVolume)
{
    // (5.45) Starting Time
    float t0 = Mathf.Pow((2 * roomVolume * Mathf.Log(2)) / (4 * Mathf.PI * Mathf.Pow(speedOfSound, 3)), 1.0f / 3.0f);
    float t = t0;

    List<float> diracPulsesPoissonProcess = new List<float>();
    List<float> timeStepValues = new List<float>();
    int i = 0;
    // Sequence of dirac deltas, mimic density of reflection in an IR of certain volume, modelled as a temporal Poisson process 
    while (t < impulseResponseLength)
    {
        if (i >= 100000) // SAFETY FIRST
            break;
        i++;

        timeStepValues.Add(t);

        // Dirac deltas falling on the latter half of a sampling interval are taken to be negative-valued
        if (Mathf.Round(t * samplingFrequency) - (t * samplingFrequency) < 0)
            diracPulsesPoissonProcess.Add(1);
        else
            diracPulsesPoissonProcess.Add(-1);

        // (5.44) mean event occurrence µ  
        float meanEventOccurence = (4 * Mathf.PI * speedOfSound * speedOfSound * speedOfSound * t * t) / roomVolume;
        if (meanEventOccurence >= 10000) // 10.000 hz is enough resolution
            meanEventOccurence = 10000;

        // Debug.Log("t: " + t + " µ: " + mu + " i: " + i);

        float z = UnityEngine.Random.Range(0.00001f, 1.0f); //  0 < z <= 1
        // (5.43)  time between one pulse and another
        float deltaTevent = Mathf.Log(1 / z) / meanEventOccurence;
        t = t + deltaTevent;
    }

    // Debug.Log(timeStepValues.Count + " " + diracPulsesPoissonProcess.Count); // Both have length of ~< 10000

    // Copy the dirac pulses to fit into impulseresponse length * samplingFrequency - Is this oversampling?
    float[] diracSequenceIRLength = new float[(int)Mathf.Ceil(impulseResponseLength * samplingFrequency)]; // 44100 samples

    // RandomSequence is a sequence of dirac pulses, where the dirac pulses are placed at the timeStepValues
    for (int j = 0; j < timeStepValues.Count; j++)
    {
        int index = Mathf.FloorToInt(timeStepValues[j] * samplingFrequency); // 0...44100
        diracSequenceIRLength[index] = diracPulsesPoissonProcess[j]; // this may override values if index is the same (but that is ok?!) 
    }

    // WriteFileFromArray.WriteCSV(diracSequenceIRLength, "1_diracSequenceIRLength" + roomVolume);

    // Why 8192? For better resolution in frequency doamin. Number of Spectral Lines = FFT Size / 2 + 1
    //The FFT Size  defines the number of frequency bins in the spectrogram.
    // Range (Nyquist): 0 - SamplingRate / 2 
    int nFFT = 8192; // FFT Size defines the frequency resolution of the spectrogram. 

    // 882 samples correspond to a window duration of approximately 20 milliseconds at a sampling frequency of 44.1 kHz
    int hannWindowLength = 882; // The window length defines the time resolution of the spectrogram. 
    // Also the window length defines the number of time frames in the spectrogram.
    // frames = (sampleCount - windowLength) / hopSize + 1

    // Hop size for 50% overlap
    int hopSize = hannWindowLength / 2; // The Hopsize defines the overlap between two consecutive windows. The hopsize is usually chosen to be 50% of the window length.

    float frequency_resolution = (float)samplingFrequency / nFFT;
    int frequency_bins = (nFFT / 2) + 1;

    // 2. and filtered by an asymmetrical high/low-pass-combination with differing slope and shape of a Raised Cosine Filter RCF (n) with
    // Raised Cosine Filter RCF (5.46)
    float[, ] RCF = new float[frequencyBands, frequency_bins]; // 6 x 4097 = 24582 // Filter für die einzelnen Frequenzbänder pro Histogram

    // RCF is correctly (checked against MATLAB)
    for (int iFreqBand = 0; iFreqBand < frequencyBands; iFreqBand++)
    {
        for (int iFreq_bin = 0; iFreq_bin < frequency_bins; iFreq_bin++)
        {
            float frequency = iFreq_bin * frequency_resolution;

            if (frequency < frequencyCenter[iFreqBand] && frequency >= lowerCutoffFrequency[iFreqBand])
                RCF[iFreqBand, iFreq_bin] = 0.5f * (1 + Mathf.Cos(2 * Mathf.PI * frequency / frequencyCenter[iFreqBand]));
            else if (frequency < upperCutoffFrequency[iFreqBand] && frequency >= frequencyCenter[iFreqBand])
                RCF[iFreqBand, iFreq_bin] = 0.5f * (1 - Mathf.Cos(2 * Mathf.PI * frequency / (frequencyCenter[iFreqBand] + 1)));
        }
    }

    // Frequency Bins = FFT Size / 2 + 1 
    Stft stft = new Stft(hannWindowLength, hopSize, NWaves.Windows.WindowType.Hann, nFFT);
    // List<float[]> spectograms = stft.Spectrogram(diracPulsesPoissonProcess.ToArray());
    // Debug.Log(spectograms.Count); // Why 22 Specotgrams? Because 22 Frames 

    // Filter the Poisson Sequence with the RCF - 6 Bandpassfilters
    int frameLength = hopSize; // 441
    int numberOfFrames = diracSequenceIRLength.Length / frameLength; // 44100 / 441 = 100 frames
    float[, ] filteredDiracSequence = new float[diracSequenceIRLength.Length, frequencyBands]; // 44100 x 6

    // Iterate over the diracSequence via frames and each frame is filtered by the RCF
    for (int frame = 1; frame <= numberOfFrames; frame++)
    {
        int from = (frame - 1) * frameLength + 1;
        int to = frame * frameLength;

        float[] binsPerFrame = new float[frameLength];
        // Debug.Log("Frame " + frame + " from: " + from + " to: " + to + " " + (to - from));

        // Pick Frame from PulseSequence: Frame 1 from: 0 to: 441 - Frame 100 from: 43659 to: 44100 Length: 441
        for (int sequence = from; sequence < to; sequence++)
            binsPerFrame[sequence - from] = diracSequenceIRLength[sequence];

        // Frequency spectrum of the frame of the Pulse Sequence
        // 1. Consequently, the Dirac delta sequence hd (see Fig. 5.18(a)) is transformed to frequency domain (see Fig. 5.18(b)) 
        // Transform to frequency domain - short-time Fourier transform and inverse short-time Fourier transform (STFT/ISTFT)
        List < (float[], float[]) > stftFrameResult = stft.Direct(binsPerFrame); // Lengt 8192 - Values Index: 0...4096, Access until 4097. Similiar to Matlab

        // Filter with RCF
        List < (float[], float[]) > resultElementWise = ElementWiseMultiplicationRCF(stftFrameResult, RCF); // 6 x (2x) 4097 - Similiar to Matlab

        // Create 6 lists for each frequency band for the inverse STFT
        List < (float[], float[]) > stftFrameResultBand1 = new List < (float[], float[]) > ();
        List < (float[], float[]) > stftFrameResultBand2 = new List < (float[], float[]) > ();
        List < (float[], float[]) > stftFrameResultBand3 = new List < (float[], float[]) > ();
        List < (float[], float[]) > stftFrameResultBand4 = new List < (float[], float[]) > ();
        List < (float[], float[]) > stftFrameResultBand5 = new List < (float[], float[]) > ();
        List < (float[], float[]) > stftFrameResultBand6 = new List < (float[], float[]) > ();

        // Add the result of the elementwise multiplication to the lists for each FrequencyBand
        stftFrameResultBand1.Add(resultElementWise[0]); // Length 1
        stftFrameResultBand2.Add(resultElementWise[1]);
        stftFrameResultBand3.Add(resultElementWise[2]);
        stftFrameResultBand4.Add(resultElementWise[3]);
        stftFrameResultBand5.Add(resultElementWise[4]);
        stftFrameResultBand6.Add(resultElementWise[5]);

        // Inverse STFT for each FrequencyBand
        float[] inverseBand1 = stft.Inverse(stftFrameResultBand1, false); // 8633 Array Size: Values to Index 880
        float[] inverseBand2 = stft.Inverse(stftFrameResultBand2, false);
        float[] inverseBand3 = stft.Inverse(stftFrameResultBand3, false);
        float[] inverseBand4 = stft.Inverse(stftFrameResultBand4, false);
        float[] inverseBand5 = stft.Inverse(stftFrameResultBand5, false);
        float[] inverseBand6 = stft.Inverse(stftFrameResultBand6, false);

        // Then the filtered spectrum is transformed back to the time domain resulting in n bandpass-filtered Dirac delta sequences hd(n)
        int inverseSamplePoint = 0;
        for (int samplePoint = from; samplePoint < to; samplePoint++) // 441 or 881?! Matlab example uses 441, nwaves returns 881 TODO: Check why why why 
        {
            filteredDiracSequence[samplePoint, 0] = inverseBand1[inverseSamplePoint]; // The filtered Sequence contains 0 values between the frames TODO: Check where the 0 are coming from.
            filteredDiracSequence[samplePoint, 1] = inverseBand2[inverseSamplePoint];
            filteredDiracSequence[samplePoint, 2] = inverseBand3[inverseSamplePoint];
            filteredDiracSequence[samplePoint, 3] = inverseBand4[inverseSamplePoint];
            filteredDiracSequence[samplePoint, 4] = inverseBand5[inverseSamplePoint];
            filteredDiracSequence[samplePoint, 5] = inverseBand6[inverseSamplePoint];
            inverseSamplePoint++;
        }
    }

    return filteredDiracSequence;
}