Как восстановить результат Сигнал-STFT-ISTFT-Сигнал, используя перекрывающиеся окна?

Извиняюсь, если ответ есть где-то на сайте (не смог найти). Я — любитель, который пытается загрузить WAV-файл, получить его амплитуду и фазовые данные (для модификации), сгенерировать спектрограмму, а затем сохранить ее как новый WAV-файл.
Я использую C ++ (Qt) и библиотеку FFTW.

Моя проблема в том, что результирующий WAV отличается от оригинала, даже если не вносить никаких изменений. Если операции FFT выполняются на всей последовательности семплов, то она выглядит точно так же, как и оригинал. Но я должен использовать STFT с перекрывающимися окнами. В этом случае я получаю искажения, приводящие к периодическим звукам трещин / удушения, также значительно изменяется форма звука.

Это можно увидеть в следующих примерах (в Audacity):

оригинал / обработан одним куском:
оригинал

обработано (windowSize = 2048, hopSize = 1024, нет оконной функции):
обработано ws = 2048, hs = 1024, wf = нет

Я не могу опубликовать больше примеров с моей репутацией, но выполнение оконной функции Хэмминга после ISTFT (а не до STFT) с помощью метода, который я использую для объединения полученных оконных сэмплов, дает хороший звук. Но форма волны все еще довольно различна, в основном наблюдаются значительные потери в пиках.

Я думаю, проблема в том, как я комбинирую результат ISTFT в новую последовательность образцов. Как правильно это сделать? Пример в C ++ был бы очень признателен.

• РЕДАКТИРОВАТЬ *
  Как правильно указано SleuthEye Я сделал ошибку в коде.
  Код скорректирован. Форма волны и звук кажутся идеальными сейчас даже без применения оконной функции. Тем не менее, метод правильный для такой операции?

Вот соответствующий источник:

// getSampleNormalized(n) returns sample n of 1 channel in -1.0 to 1.0 range
// getSampleCount()      returns sample count of 1 channel
// quint32 is just unsigned int

quint32 windowSize     = 2048;
quint32 windowSizeHalf = windowSize / 2 + 1;
quint32 slideWindowBy  = 1024; // hopSize
quint32 windowCount    = getSampleCount() / slideWindowBy;
if ( (windowCount * slideWindowBy) < getSampleCount()){
windowCount += 1;
}
quint32 newSampleCount = windowCount * slideWindowBy + ( windowSize - slideWindowBy );

double *window          = new double[windowSize];
fftw_complex *fftResult = new fftw_complex[windowSizeHalf];
fftw_complex *fftWindow = new fftw_complex[windowSizeHalf];
double *result          = new double[windowSize];

double **magnitudes    = new double*[windowCount];
double **phases        = new double*[windowCount];
double **signalWindows = new double*[windowCount];
for (int i = 0; i < windowCount; ++i){
magnitudes[i]    = new double[windowSizeHalf];
phases[i]        = new double[windowSizeHalf];
signalWindows[i] = new double[windowSize];
}

double *sampleSignals = new double[newSampleCount];

fftw_plan fftPlan  = fftw_plan_dft_r2c_1d( windowSize, window,    fftResult, FFTW_ESTIMATE );
fftw_plan ifftPlan = fftw_plan_dft_c2r_1d( windowSize, fftWindow, result,    FFTW_ESTIMATE );

// STFT
for ( int currentWindow = 0; currentWindow < windowCount; ++currentWindow ){
for (int i = 0; i < windowSize; ++i){
quint32 currentSample = currentWindow * slideWindowBy + i;
if ( ( currentSample ) < getSampleCount() ){
window[i] = getSampleNormalized( currentSample ); // * ( windowHamming( i, windowSize ) );
}
else{
window[i] = 0.0;
}
}

fftw_execute(fftPlan);

for (int i = 0; i < windowSizeHalf; ++i){
magnitudes[currentWindow][i] = sqrt( fftResult[i][0]*fftResult[i][0] + fftResult[i][1]*fftResult[i][1] );
phases[currentWindow][i]     = atan2( fftResult[i][1], fftResult[i][0] );
}
}

// INVERSE STFT
for ( int currentWindow = 0; currentWindow < windowCount; ++currentWindow ){
for ( int i = 0; i < windowSizeHalf; ++i ){
fftWindow[i][0] = magnitudes[currentWindow][i] * cos( phases[currentWindow][i] );  // Real
fftWindow[i][1] = magnitudes[currentWindow][i] * sin( phases[currentWindow][i] );  // Imaginary
}

fftw_execute(ifftPlan);

for ( int i = 0; i < windowSize; ++i ){
signalWindows[currentWindow][i] = result[i] / windowSize;            // getting normalized result
//signalWindows[currentWindow][i] *= (windowHamming( i, windowSize )); // applying Hamming window function
}
}

quint32 pos;

// HERE WE COMBINE RESULTED WINDOWS

// COMBINE AND AVERAGE
// 1st window should be full replace
for ( int i = 0; i < windowSize; ++i ){
sampleSignals[i] = signalWindows[0][i];
}
// 2nd window and onwards: combine with previous ones
for ( int currentWindow = 1; currentWindow < windowCount; ++currentWindow ){
// combine and average with data from previous window
for ( int i = 0; i < (windowSize - slideWindowBy); ++i ){
pos = currentWindow * slideWindowBy + i;
sampleSignals[pos] = (sampleSignals[pos] + signalWindows[currentWindow][i]) * 0.5;
}
// simply replace for the rest
for ( int i = (windowSize - slideWindowBy); i < windowSize; ++i ){
pos = currentWindow * slideWindowBy + i;
sampleSignals[pos] = signalWindows[currentWindow][i];
}
}

// then just save the wav file...