|
| 1 | +#include "dsp.hpp" |
| 2 | +#include <algorithm> |
| 3 | +#include <cmath> |
| 4 | +#include <cstdlib> |
| 5 | +#include <iostream> |
| 6 | +#include <memory> |
| 7 | +#include <string> |
| 8 | +#include <unsupported/Eigen/FFT> |
| 9 | +#include <vector> |
| 10 | + |
| 11 | +// forward declaration of inner stft |
| 12 | +void stft_inner(struct demucsonnx::stft_buffers &stft_buf, |
| 13 | + Eigen::FFT<float> &cfg); |
| 14 | + |
| 15 | +void istft_inner(struct demucsonnx::stft_buffers &stft_buf, |
| 16 | + Eigen::FFT<float> &cfg); |
| 17 | + |
| 18 | +// reflect padding |
| 19 | +void pad_signal(struct demucsonnx::stft_buffers &stft_buf) |
| 20 | +{ |
| 21 | + // copy from stft_buf.padded_waveform_mono_in+pad into stft_buf.pad_start, |
| 22 | + // stft_buf.pad_end |
| 23 | + std::copy_n(stft_buf.padded_waveform_mono_in.begin() + stft_buf.pad, |
| 24 | + stft_buf.pad, stft_buf.pad_start.begin()); |
| 25 | + std::copy_n(stft_buf.padded_waveform_mono_in.end() - 2 * stft_buf.pad, |
| 26 | + stft_buf.pad, stft_buf.pad_end.begin()); |
| 27 | + |
| 28 | + std::reverse(stft_buf.pad_start.begin(), stft_buf.pad_start.end()); |
| 29 | + std::reverse(stft_buf.pad_end.begin(), stft_buf.pad_end.end()); |
| 30 | + |
| 31 | + // copy stft_buf.pad_start into stft_buf.padded_waveform_mono_in |
| 32 | + std::copy_n(stft_buf.pad_start.begin(), stft_buf.pad, |
| 33 | + stft_buf.padded_waveform_mono_in.begin()); |
| 34 | + |
| 35 | + // copy stft_buf.pad_end into stft_buf.padded_waveform_mono_in |
| 36 | + std::copy_n(stft_buf.pad_end.begin(), stft_buf.pad, |
| 37 | + stft_buf.padded_waveform_mono_in.end() - stft_buf.pad); |
| 38 | +} |
| 39 | + |
| 40 | +Eigen::FFT<float> get_fft_cfg() |
| 41 | +{ |
| 42 | + Eigen::FFT<float> cfg; |
| 43 | + |
| 44 | + cfg.SetFlag(Eigen::FFT<float>::Speedy); |
| 45 | + // cfg.SetFlag(Eigen::FFT<float>::HalfSpectrum); |
| 46 | + // cfg.SetFlag(Eigen::FFT<float>::Unscaled); |
| 47 | + |
| 48 | + return cfg; |
| 49 | +} |
| 50 | + |
| 51 | +void demucsonnx::stft( |
| 52 | + struct stft_buffers &stft_buf, |
| 53 | + const Eigen::MatrixXf &waveform, |
| 54 | + Eigen::Tensor3dXcf &spec) |
| 55 | +{ |
| 56 | + // get the fft config |
| 57 | + Eigen::FFT<float> cfg = get_fft_cfg(); |
| 58 | + |
| 59 | + /*****************************************/ |
| 60 | + /* operate on each channel sequentially */ |
| 61 | + /*****************************************/ |
| 62 | + |
| 63 | + for (int channel = 0; channel < 2; ++channel) |
| 64 | + { |
| 65 | + Eigen::VectorXf row_vec = waveform.row(channel); |
| 66 | + |
| 67 | + std::copy_n(row_vec.data(), row_vec.size(), |
| 68 | + stft_buf.padded_waveform_mono_in.begin() + stft_buf.pad); |
| 69 | + |
| 70 | + // apply padding equivalent to center padding with center=True |
| 71 | + // in torch.stft: |
| 72 | + // https://pytorch.org/docs/stable/generated/torch.stft.html |
| 73 | + |
| 74 | + // reflect pads stft_buf.padded_waveform_mono in-place |
| 75 | + pad_signal(stft_buf); |
| 76 | + |
| 77 | + // does forward fft on stft_buf.padded_waveform_mono, stores spectrum in |
| 78 | + // complex_spec_mono |
| 79 | + stft_inner(stft_buf, cfg); |
| 80 | + |
| 81 | + for (int i = 0; i < stft_buf.nb_bins; ++i) |
| 82 | + { |
| 83 | + for (int j = 0; j < stft_buf.nb_frames; ++j) |
| 84 | + { |
| 85 | + spec(channel, i, j) = stft_buf.complex_spec_mono[j][i]; |
| 86 | + } |
| 87 | + } |
| 88 | + } |
| 89 | +} |
| 90 | + |
| 91 | +void demucsonnx::istft( |
| 92 | + struct stft_buffers &stft_buf, |
| 93 | + const Eigen::Tensor3dXcf &spec, |
| 94 | + Eigen::MatrixXf &waveform) |
| 95 | +{ |
| 96 | + // get the fft config |
| 97 | + Eigen::FFT<float> cfg = get_fft_cfg(); |
| 98 | + |
| 99 | + /*****************************************/ |
| 100 | + /* operate on each channel sequentially */ |
| 101 | + /*****************************************/ |
| 102 | + |
| 103 | + for (int channel = 0; channel < 2; ++channel) |
| 104 | + { |
| 105 | + // Populate the nested vectors |
| 106 | + for (int i = 0; i < stft_buf.nb_bins; ++i) |
| 107 | + { |
| 108 | + for (int j = 0; j < stft_buf.nb_frames; ++j) |
| 109 | + { |
| 110 | + stft_buf.complex_spec_mono[j][i] = spec(channel, i, j); |
| 111 | + } |
| 112 | + } |
| 113 | + |
| 114 | + // does inverse fft on stft_buf.complex_spec_mono, stores waveform in |
| 115 | + // padded_waveform_mono |
| 116 | + istft_inner(stft_buf, cfg); |
| 117 | + |
| 118 | + // copies waveform_mono into stft_buf.waveform past first pad samples |
| 119 | + waveform.row(channel) = Eigen::Map<Eigen::MatrixXf>( |
| 120 | + stft_buf.padded_waveform_mono_out.data() + stft_buf.pad, 1, |
| 121 | + stft_buf.padded_waveform_mono_out.size() - FFT_WINDOW_SIZE); |
| 122 | + } |
| 123 | +} |
| 124 | + |
| 125 | +void stft_inner(struct demucsonnx::stft_buffers &stft_buf, |
| 126 | + Eigen::FFT<float> &cfg) |
| 127 | +{ |
| 128 | + int frame_idx = 0; |
| 129 | + |
| 130 | + // Loop over the waveform with a stride of hop_size |
| 131 | + for (std::size_t start = 0; |
| 132 | + start <= |
| 133 | + stft_buf.padded_waveform_mono_in.size() - demucsonnx::FFT_WINDOW_SIZE; |
| 134 | + start += demucsonnx::FFT_HOP_SIZE) |
| 135 | + { |
| 136 | + // Apply window and run FFT |
| 137 | + for (int i = 0; i < demucsonnx::FFT_WINDOW_SIZE; ++i) |
| 138 | + { |
| 139 | + stft_buf.windowed_waveform_mono[i] = |
| 140 | + stft_buf.padded_waveform_mono_in[start + i] * |
| 141 | + stft_buf.window[i]; |
| 142 | + } |
| 143 | + cfg.fwd(stft_buf.complex_spec_mono[frame_idx], |
| 144 | + stft_buf.windowed_waveform_mono); |
| 145 | + // now scale stft_buf.complex_spec_mono[frame_idx] by 1.0f / |
| 146 | + // sqrt(float(FFT_WINDOW_SIZE))) |
| 147 | + |
| 148 | + for (int i = 0; i < demucsonnx::FFT_WINDOW_SIZE / 2 + 1; ++i) |
| 149 | + { |
| 150 | + stft_buf.complex_spec_mono[frame_idx][i] *= |
| 151 | + 1.0f / sqrt(float(demucsonnx::FFT_WINDOW_SIZE)); |
| 152 | + } |
| 153 | + frame_idx++; |
| 154 | + } |
| 155 | +} |
| 156 | + |
| 157 | +void istft_inner(struct demucsonnx::stft_buffers &stft_buf, |
| 158 | + Eigen::FFT<float> &cfg) |
| 159 | +{ |
| 160 | + // clear padded_waveform_mono |
| 161 | + std::fill(stft_buf.padded_waveform_mono_out.begin(), |
| 162 | + stft_buf.padded_waveform_mono_out.end(), 0.0f); |
| 163 | + |
| 164 | + // Loop over the input with a stride of (hop_size) |
| 165 | + for (int start = 0; start < stft_buf.nb_frames * demucsonnx::FFT_HOP_SIZE; |
| 166 | + start += demucsonnx::FFT_HOP_SIZE) |
| 167 | + { |
| 168 | + int frame_idx = start / demucsonnx::FFT_HOP_SIZE; |
| 169 | + // undo sqrt(nfft) scaling |
| 170 | + for (int i = 0; i < demucsonnx::FFT_WINDOW_SIZE / 2 + 1; ++i) |
| 171 | + { |
| 172 | + stft_buf.complex_spec_mono[frame_idx][i] *= |
| 173 | + sqrt(float(demucsonnx::FFT_WINDOW_SIZE)); |
| 174 | + } |
| 175 | + // Run iFFT |
| 176 | + cfg.inv(stft_buf.windowed_waveform_mono, |
| 177 | + stft_buf.complex_spec_mono[frame_idx]); |
| 178 | + |
| 179 | + // Apply window and add to output |
| 180 | + for (int i = 0; i < demucsonnx::FFT_WINDOW_SIZE; ++i) |
| 181 | + { |
| 182 | + // x[start+i] is the sum of squared window values |
| 183 | + // https://github.com/librosa/librosa/blob/main/librosa/core/spectrum.py#L613 |
| 184 | + // 1e-8f is a small number to avoid division by zero |
| 185 | + stft_buf.padded_waveform_mono_out[start + i] += |
| 186 | + stft_buf.windowed_waveform_mono[i] * stft_buf.window[i] * 1.0f / |
| 187 | + float(demucsonnx::FFT_WINDOW_SIZE) / |
| 188 | + (stft_buf.normalized_window[start + i] + 1e-8f); |
| 189 | + } |
| 190 | + } |
| 191 | +} |
0 commit comments