音频处理中重叠相加法实现故障排查与修复请求

我是数字信号处理（DSP）新手，尝试用**重叠相加法（Overlap-Add Method）**处理音频信号，以解决调制/解调等操作后出现的信号不连续问题。但当前实现产生了意外的失真音效，信号不连续问题也未得到解决。

已尝试的方案：

• 基于资料初步实现重叠相加法，直接应用于处理后的音频信号以平滑不连续点；
• 调整操作顺序，将重叠相加步骤从信号处理后移至解调后，聚焦于合并信号实部；
• 调整缓冲区大小与重叠比例，尝试消除音频伪影。

附相关实现代码：

#include <SoapySDR/Device.h>
#include <SoapySDR/Formats.h>
#include <complex.h>
#include <fftw3.h>
#include <math.h>
#include <sndfile.h>
#include <stdio.h>   //printf
#include <stdlib.h>  //free
#include <string.h>
#define M_PI 3.14159265358979323846264338327950288
#define BUFFER_SIZE 32768
#define SAMPLE_RATE 2.4e6
#define CENTER_FREQUENCY 100e6
#define TARGET_FREQUENCY 100.3e6
#define BANDWIDTH 200e3
#define PI 3.14159265358979323846

complex float prev_sample = 0;  // Initialisieren mit 0 für den ersten Aufruf
complex float buffer[BUFFER_SIZE];
complex float output_buffer[BUFFER_SIZE];
complex float output_buffer_prev[BUFFER_SIZE];
float demodulated_sample[BUFFER_SIZE];
float demodulated_sample_prev[BUFFER_SIZE];
size_t buffer_index = 0;

void add_complex_arr(complex float *arr1, complex float *arr2, size_t len) {
  for (size_t i = 0; i < len; i++) {
    arr1[i] += arr2[i];
  }
}

void generate_hann_window(float *window, size_t length) {
  for (size_t i = 0; i < length; ++i) {
    window[i] = 0.5 * (1 - cos(2 * M_PI * i / (length - 1)));
  }
}


// Global variables for FFTW plans and arrays
fftw_complex *in_global, *out_global;
fftw_plan p_global, p_inv_global;

void initialize_fftw(size_t length) {
  // Allocate memory for FFTW input and output arrays
  in_global = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * length);
  out_global = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * length);

  // Create FFTW plans
  p_global = fftw_plan_dft_1d(length, in_global, out_global, FFTW_FORWARD,
                              FFTW_ESTIMATE);
  p_inv_global = fftw_plan_dft_1d(length, out_global, in_global, FFTW_BACKWARD,
                                  FFTW_ESTIMATE);
}

void cleanup_fftw() {
  // Cleanup
  fftw_destroy_plan(p_global);
  fftw_destroy_plan(p_inv_global);
  fftw_free(in_global);
  fftw_free(out_global);
}

// Big credits to HA7ILM for the FM Demodulation - This is very effecient and
// sounds great! https://github.com/ha7ilm/csdr
#define fmdemod_quadri_K 0.340447550238101026565118445432744920253753662109375
// this constant ensures proper scaling for qa_fmemod testcases for SNR
// calculation and more.
void fm_demodulate(complex float *in, float *out, size_t len,
                   complex float last_sample) {
  if (len == 10) return;

  float temp_dq[len];
  float temp_di[len];
  float temp[len];

  // Handle dq
  temp_dq[0] = cimagf(in[0]) - cimagf(last_sample);
  for (size_t i = 1; i < len; i++) {
    temp_dq[i] = cimagf(in[i]) - cimagf(in[i - 1]);
  }

  // Handle di
  temp_di[0] = crealf(in[0]) - crealf(last_sample);
  for (size_t i = 1; i < len; i++) {
    temp_di[i] = crealf(in[i]) - crealf(in[i - 1]);
  }

  // Output numerator
  for (size_t i = 0; i < len; i++) {
    out[i] = crealf(in[i]) * temp_dq[i] - cimagf(in[i]) * temp_di[i];
  }

  // Output denominator
  for (size_t i = 0; i < len; i++) {
    temp[i] = crealf(in[i]) * crealf(in[i]) + cimagf(in[i]) * cimagf(in[i]);
  }

  // Output division
  for (size_t i = 0; i < len; i++) {
    out[i] = (temp[i]) ? fmdemod_quadri_K * out[i] / temp[i] : 0;
  }
}

int process_block(fftw_complex *fft_result, size_t length,
                  complex float *output, float target_frequency) {
  // Calculate indices for the target frequency band
  float bin_width = SAMPLE_RATE / length;
  size_t start_bin = (size_t)((target_frequency - CENTER_FREQUENCY -
                               BANDWIDTH / 2 + SAMPLE_RATE) /
                              bin_width) %
                     length;
  size_t end_bin = (size_t)((target_frequency - CENTER_FREQUENCY +
                             BANDWIDTH / 2 + SAMPLE_RATE) /
                            bin_width) %
                   length;

  // Allocate memory for the reduced FFT results containing only the relevant
  // bins
  size_t relevant_bins = end_bin - start_bin + 1;
  fftw_complex *relevant_fft =
      (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * relevant_bins);

  // Copy only the relevant FFT results
  for (size_t i = start_bin, j = 0; i <= end_bin; i++, j++) {
    relevant_fft[j][0] = fft_result[i][0];
    relevant_fft[j][1] = fft_result[i][1];
  }

  // Create a temporary plan for inverse FFT on the reduced set of bins
  fftw_plan p_relevant_inv = fftw_plan_dft_1d(
      relevant_bins, relevant_fft, relevant_fft, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_execute(p_relevant_inv);  // Execute the inverse FFT

  // Normalize and copy the IFFT result to the output buffer
  for (size_t i = 0; i < relevant_bins; i++) {
    output[i] = relevant_fft[i][0] / relevant_bins +
                relevant_fft[i][1] * I / relevant_bins;
  }

  add_complex_arr((complex float *)output, output_buffer_prev,
                  relevant_bins / 2);

  size_t half_index = relevant_bins / 2;
  // Copy the second half of demodulated_sample to demodulated_sample_prev
  memcpy(output_buffer_prev, output + half_index, half_index * sizeof(float));

  // Cleanup
  fftw_destroy_plan(p_relevant_inv);
  fftw_free(relevant_fft);

  return relevant_bins;  // Return the number of samples in the output buffer
}

void execute_fft(complex float *input, size_t length, fftw_complex *output) {
  // Copy input data to FFTW input array
  for (size_t i = 0; i < length; i++) {
    in_global[i][0] = crealf(input[i]);
    in_global[i][1] = cimagf(input[i]);
  }

  // Execute FFT
  fftw_execute(p_global);

  // Copy FFT result to output array
  memcpy(output, out_global, sizeof(fftw_complex) * length);
}

void fm_demodulate_real(const float *in, float *out, size_t len) {
  // Ensure there's at least two samples to process
  if (len < 2) return;

  // Iterate through each sample, except the last one
  for (size_t i = 0; i < len - 1; ++i) {
    // Simple difference equation as an approximation
    // This is a basic form of differentiation; you may need to adjust scaling
    out[i] = in[i + 1] - in[i];
  }

  // Handle the last sample manually or replicate the previous difference
  out[len - 1] = out[len - 2];
}

// Main processing loop
void process_data(complex float *data, size_t total_size, SNDFILE *file) {
  complex float windowed_data[total_size];
  complex float output_data[total_size];
  fftw_complex *fft_result =
      (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * BUFFER_SIZE);

  float hann_window[total_size];
  generate_hann_window(hann_window, total_size);
  // Apply the Hann window to the output_buffer
  for (size_t i = 0; i < total_size; i++) {
    data[i] *= hann_window[i];
  }

  execute_fft(data, total_size, fft_result);
  int received_data =
      process_block(fft_result, total_size, output_buffer, TARGET_FREQUENCY);

  // Calculate the start index for the second half of the demodulated_sample

  fm_demodulate(output_buffer, demodulated_sample, received_data, prev_sample);

  prev_sample = output_buffer[received_data / 2 - 1];

  sf_writef_float(file, demodulated_sample, received_data);
}

int main(void) {
  size_t length;

  // enumerate devices
  SoapySDRKwargs *results = SoapySDRDevice_enumerate(NULL, &length);
  for (size_t i = 0; i < length; i++) {
    printf("Found device #%d: ", (int)i);
    for (size_t j = 0; j < results[i].size; j++) {
      printf("%s=%s, ", results[i].keys[j], results[i].vals[j]);
    }
    printf("\n");
  }
  SoapySDRKwargsList_clear(results, length);

  // create device instance
  SoapySDRKwargs args = {};
  SoapySDRKwargs_set(&args, "driver", "rtlsdr");

  SoapySDRDevice *sdr = SoapySDRDevice_make(&args);
  SoapySDRKwargs_clear(&args);

  if (sdr == NULL) {
    printf("SoapySDRDevice_make fail: %s\n", SoapySDRDevice_lastError());
    return EXIT_FAILURE;
  }

  // apply settings
  if (SoapySDRDevice_setSampleRate(sdr, SOAPY_SDR_RX, 0, 2.4e6) != 0) {
    printf("setSampleRate fail: %s\n", SoapySDRDevice_lastError());
  }
  if (SoapySDRDevice_setFrequency(sdr, SOAPY_SDR_RX, 0, 100e6, NULL) != 0) {
    printf("setFrequency fail: %s\n", SoapySDRDevice_lastError());
  }

  // setup a stream (complex floats)
  SoapySDRStream *rxStream = SoapySDRDevice_setupStream(
      sdr, SOAPY_SDR_RX, SOAPY_SDR_CF32, NULL, 0, NULL);
  if (rxStream == NULL) {
    printf("setupStream fail: %s\n", SoapySDRDevice_lastError());
    SoapySDRDevice_unmake(sdr);
    return EXIT_FAILURE;
  }
  SoapySDRDevice_activateStream(sdr, rxStream, 0, 0, 0);  // start streaming

  // create a re-usable buffer for rx samples
  complex float buff[32768];

  SF_INFO info;
  info.samplerate = 200e3;
  info.channels = 1;
  info.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
  SNDFILE *file = sf_open("output.wav", SFM_WRITE, &info);

  initialize_fftw(BUFFER_SIZE);

  // receive some samples
  for (size_t i = 0; i < 500; i++) {
    void *buffs[] = {buff};  // array of buffers
    int flags;               // flags set by receive operation
    long long timeNs;        // timestamp for receive buffer
    int ret = SoapySDRDevice_readStream(sdr, rxStream, buffs, 32768, &flags,
                                        &timeNs, 1000000);

    // Check for errors in readStream
    if (ret < 0) {
      fprintf(stderr, "Error reading stream: %s\n", SoapySDR_errToStr(ret));
      continue;  // Skip to next iteration on error
    }
    // Copy received samples from buff to buffer
    // memcpy(buffer + buffer_index, buff, ret * sizeof(complex float));
    // buffer_index += ret;
    // if (buffer_index >= BUFFER_SIZE) {

    process_data(buff, BUFFER_SIZE, file);

    // buffer_index = 0;
    //}
  }

  sf_close(file);

  // shutdown the stream
  SoapySDRDevice_deactivateStream(sdr, rxStream, 0, 0);  // stop streaming
  SoapySDRDevice_closeStream(sdr, rxStream);

  // cleanup device handle
  SoapySDRDevice_unmake(sdr);

  printf("Done\n");
  return EXIT_SUCCESS;
}

内容的提问来源于stack exchange，提问作者StefenMA