[C++ 基于Eigen库实现CRN前向推理]
第一部分:WavFile.class (实现读取wav/pcm,实现STFT)
- 前言:(Eigen库使用记录)
- 第一部分:WavFile.class (实现读取wav/pcm,实现STFT)
- 第二部分:Conv2d实现
- 第三部分:TransposedConv2d实现 (mimo,padding,stride,dilation,kernel,outpadding)
- 第四部分:NonLinearity (Sigmoid,Tanh,ReLU,ELU,Softplus)
- 第五部分:LSTM
- GITHUB仓库
1. 读写语音文件(wav/pcm)
1.1 wav文件读写
语音文件有很多后缀格式,这里只实现了两种后缀名,.wav和.pcm。
- WAV文件遵守资源交换文件格式之规则,在文件的前44(或46)字节放置标头(header),使播放器或编辑器能够简单掌握文件的基本信息,其内容以区块(chunk)为最小单位,每一区块长度为4字节,而区块之上则由子区块包裹,每一子区块长度不拘,但须在前头先宣告标签及长度(字节)。
- 标头的前3个区块记录文件格式及长度;
- 接着第一个子区块包含8个区块,记录声道数量、采样率等信息;
- 接着第二个子区块才是真正的音频资料,长度则视音频长度而定。
- 内容如下表所示。须注意的是,每个区块的端序不尽相同,而音频内容本身则是采用小端序。
-
读Wav文件
图中所示是一个wav的标准格式,但wav还有其他不同的存储格式,有些有扩展块之类的,需要根据具体的语音文件来确定。比如我用python中soudfile生成的语音,其format_tag=3,format_length=18(通常为16)。由于我没用到其他格式的情况,直接在下方用if else讨论format_tag=1和=3的情况。
总而言之,读取wav其实就是一个解析头文件的过程,把对应的头信息存下来,再去读取数据就好了。-
读wav文件的具体代码给出(后面会给整体的h文件和cpp文件):void Wav_File::LoadWavFile(const char *file_path) { FILE *fp = fopen(file_path, "rb"); if (fp) { fread(this->id_riff, sizeof(char), 4, fp); this->id_riff[4] = '\0'; fread(&this->file_size, sizeof(int16_t), 2, fp); fread(this->id_wave, sizeof(char), 4, fp); this->id_wave[4] = '\0'; fread(this->id_fmt, sizeof(char), 4, fp); this->id_fmt[4] = '\0'; fread(&this->format_length, sizeof(int16_t), 2, fp); fread(&this->format_tag, sizeof(int16_t), 1, fp); fread(&this->channels, sizeof(int16_t), 1, fp); fread(&this->sample_rate, sizeof(int16_t), 2, fp); fread(&this->avg_bytes_sec, sizeof(int16_t), 2, fp); fread(&this->block_align, sizeof(int16_t), 1, fp); fread(&this->bits_per_sample, sizeof(int16_t), 1, fp); if (this->format_tag == 1) { fread(this->id_data, sizeof(char), 4, fp); this->id_data[4] = '\0'; fread(&this->data_size, sizeof(int16_t), 2, fp); this->wav_size = this->data_size / sizeof(int16_t); this->ip_raw = (int16_t *) malloc(this->data_size); fread(this->ip_raw, sizeof(int16_t), this->wav_size, fp); this->data = (float_t *) malloc(this->wav_size * sizeof(float_t)); for (int i = 0; i < this->wav_size; i++) { this->data[i] = (float_t) (this->ip_raw[i] * 1.0 / 32768); } } else if (this->format_tag == 3) { fread(&this->cb_size, sizeof(int16_t), 1, fp); fread(this->id_data, sizeof(char), 4, fp); this->id_data[4] = '\0'; fread(&this->data_size, sizeof(int16_t), 2, fp); this->wav_size = this->data_size / sizeof(float); this->fp_raw = (float *) malloc(this->data_size); fread(this->fp_raw, sizeof(float), this->wav_size, fp); this->data = (float_t *) malloc(this->wav_size * sizeof(float_t)); for (int i = 0; i < this->wav_size; i++) { this->data[i] = this->fp_raw[i]; } } } fclose(fp); }
-
写出Wav文件
Wav文件写出和读入是完全对应的,把读到的信息反向写出就可以了。
如果对数据有所操作,需要去找对应的变量进行修改,比如长度,需要更新总文件的file_size以及data chunk的data_size。
具体代码给出:void Wav_File::WriteWavFile(const char *dest_path) { for (int i = 0; i < wav_size; i++) { if (format_tag == 1) this->ip_raw[i] = (int16_t) (this->data[i] * 32768); else this->fp_raw[i] = this->data[i]; } FILE *fp = fopen(dest_path, "wb"); if (fp) { fwrite(this->id_riff, sizeof(char), 4, fp); fwrite(&this->file_size, sizeof(int16_t), 2, fp); fwrite(this->id_wave, sizeof(char), 4, fp); fwrite(this->id_fmt, sizeof(char), 4, fp); fwrite(&this->format_length, sizeof(int16_t), 2, fp); fwrite(&this->format_tag, sizeof(int16_t), 1, fp); fwrite(&this->channels, sizeof(int16_t), 1, fp); fwrite(&this->sample_rate, sizeof(int16_t), 2, fp); fwrite(&this->avg_bytes_sec, sizeof(int16_t), 2, fp); fwrite(&this->block_align, sizeof(int16_t), 1, fp); fwrite(&this->bits_per_sample, sizeof(int16_t), 1, fp); if (this->format_tag == 3) { fwrite(&this->cb_size, sizeof(int16_t), 1, fp); } fwrite(this->id_data, sizeof(char), 4, fp); fwrite(&this->data_size, sizeof(int16_t), 2, fp); if (this->format_tag == 1) { fwrite(this->ip_raw, sizeof(int16_t), this->wav_size, fp); } else if (this->format_tag == 3) { fwrite(this->fp_raw, sizeof(float), this->wav_size, fp); } } fclose(fp); }
1.2 PCM文件读写
PCM文件没有格式,是直接存的数据,所以需要自己知道采样率,通道数,以及数据类型,然后直接读写就可以了,代码如下:
void Wav_File::LoadPcmFile(const char *file_path) {
FILE *fp = fopen(file_path, "rb");
// 计算文件长度
fseek(fp, 0L, SEEK_END);
this->wav_size = ftell(fp) / sizeof(float);
fseek(fp, 0L, SEEK_SET);
//申请内存空间 这边自己知道存储的是float32,所以申请float的内存空间
this->fp_raw = (float *) malloc(this->wav_size * sizeof(float));
fread(this->fp_raw, sizeof(float), this->wav_size, fp);
fclose(fp);
// 这部分是统一转换成float32存在data里,因为wav是用int16存的
for (int i = 0; i < this->wav_size; i++) {
this->data[i] = this->fp_raw[i];
}
}
void Wav_File::WritePcmFile(const char *dest_path) {
for (int i = 0; i < this->wav_size; i++) {
this->fp_raw[i] = this->data[i];
}
FILE *fp = fopen(dest_path, "wb");
fwrite(this->fp_raw, sizeof(float), this->wav_size, fp);
fclose(fp);
}
2. STFT实现
2.1 设置STFT参数,窗函数
这块主要设置窗长,窗移,窗函数类型等,并为频谱数据申请内存空间
提供了三种窗实现,hanning,hamming和sqrt窗(我自己叫的)。
可以看到有两份spec和mag,pha内存,spec是自定义复数类型,用于存储由蝶形算法生成的频谱(因为对FFT实现不是很熟悉,根据网上照猫画虎,只能对2^N窗长进行变换)。
而另一份是对应不同基的FFT,没搞懂是用的现成的,具体代码就不在这放了,可以去github上看。
void Wav_File::setSTFT(int64_t frame_size, int64_t frame_shift, const char *win_type) {
this->frame_size = frame_size;
this->frame_shift = frame_shift;
this->fft_size = this->frame_size;
this->frame_num = this->wav_size / frame_shift + 1;
this->win_type = win_type;
this->setWindow(this->win_type);
this->spec = (Complex **) malloc((this->fft_size / 2 + 1) * sizeof(*spec));
for (int i = 0; i < (this->fft_size / 2 + 1); i++)
this->spec[i] = (Complex *) malloc(this->frame_num * sizeof(Complex));
this->mag = (float_t **) malloc((this->fft_size / 2 + 1) * sizeof(*this->mag));
for (int i = 0; i < (this->fft_size / 2 + 1); i++)
this->mag[i] = (float_t *) malloc(frame_num * sizeof(float_t));
this->phase = (float_t **) malloc((this->fft_size / 2 + 1) * sizeof(*this->phase));
for (int i = 0; i < (this->fft_size / 2 + 1); i++)
this->phase[i] = (float_t *) malloc(frame_num * sizeof(float_t));
// nfft
this->spec_real = (float_t **) malloc(this->frame_num * sizeof(*spec_real));
for (int i = 0; i < this->frame_num; i++)
this->spec_real[i] = (float_t *) malloc((this->fft_size / 2 + 1) * sizeof(float_t));
this->spec_imag = (float_t **) malloc(this->frame_num * sizeof(*spec_imag));
for (int i = 0; i < this->frame_num; i++)
this->spec_imag[i] = (float_t *) malloc((this->fft_size / 2 + 1) * sizeof(float_t));
this->spec_mag = (float_t **) malloc(this->frame_num * sizeof(*spec_mag));
for (int i = 0; i < this->frame_num; i++)
this->spec_mag[i] = (float_t *) malloc((this->fft_size / 2 + 1) * sizeof(float_t));
this->spec_pha = (float_t **) malloc(this->frame_num * sizeof(*spec_pha));
for (int i = 0; i < this->frame_num; i++)
this->spec_pha[i] = (float_t *) malloc((this->fft_size / 2 + 1) * sizeof(float_t));
}
void Wav_File::setWindow(const char *winType) {
double_t a0, a1;
this->window = (float_t *) malloc(this->frame_size * sizeof(float_t));
int frac = this->frame_size - 1;
if (strcmp(winType, "hamming") == 0) {
for (int i = 0; i < this->frame_size; i++) {
this->window[i] = float(0.540000000 - 0.460000000 * cos(PI * 2 * i / frac));
}
} else if (strcmp(winType, "hanning") == 0) {
for (int i = 0; i < this->frame_size; i++) {
this->window[i] = float(0.500000000 - 0.500000000 * cos(PI * 2 * i / frac));
}
} else if (strcmp(winType, "sqrt") == 0) {
for (int i = 0; i < this->frame_size; i++) {
this->window[i] = (float) sqrtf((0.5f - 0.5f * cosf(i / (float) (this->frame_size - 1) * 2 * PI)));
}
}
}
2.2 2^N版FFT
void Wav_File::STFT() {
std::vector<float_t> padding_data;
for (int i = 0; i < this->frame_shift; i++)
padding_data.push_back(this->data[this->frame_shift - i]);
for (int i = 0; i < this->wav_size; i++)
padding_data.push_back(this->data[i]);
for (int i = 0; i < this->frame_shift; i++)
padding_data.push_back(this->data[this->wav_size - i - 2]);
this->real = (float_t *) malloc(sizeof(float_t) * this->frame_size);
this->imag = (float_t *) malloc(sizeof(float_t) * this->frame_size);
for (int i = 0; i < this->frame_num; i++) {
for (int j = 0; j < this->frame_size; j++) {
this->real[j] = padding_data[i * this->frame_shift + j] * this->window[j];
this->imag[j] = 0.0;
}
FFT(this->real, this->imag, this->fft_size);
for (int j = 0; j <= this->fft_size / 2; j++) {
this->spec[j][i].real = this->real[j];
this->spec[j][i].imag = this->imag[j];
}
}
}
void Wav_File::ISTFT() {
this->real = (float_t *) malloc(sizeof(float_t) * this->frame_size);
this->imag = (float_t *) malloc(sizeof(float_t) * this->frame_size);
auto *temp = (float_t *) malloc(sizeof(float_t) * (data_size + 2 * frame_shift));
auto *ola = (float_t *) malloc(sizeof(float_t) * (data_size + 2 * frame_shift));
for (int i = 0; i < this->frame_num; i++) {
//��N_FFT/2+1����N_FFT(����Գ�)
for (int j = 0; j < this->fft_size / 2 + 1; j++) {
real[j] = spec[j][i].real;
imag[j] = spec[j][i].imag;
}
for (int j = this->fft_size / 2 + 1; j < this->fft_size; j++) {
real[j] = spec[this->fft_size - j][i].real;
imag[j] = -spec[this->fft_size - j][i].imag;
}
IFFT(real, imag, this->fft_size);
for (int j = 0; j < frame_size; j++) {
temp[i * frame_shift + j] += (real[j] * this->window[j]);
ola[i * frame_shift + j] += (this->window[j] * this->window[j]);
}
}
for (int i = 0; i < (frame_num + 1) * frame_shift; i++) {
if (temp[i] != 0) {
temp[i] /= ola[i];
}
}
for (int i = 0; i < (frame_num - 1) * frame_shift; i++)
this->data[i] = temp[i + frame_shift];
// 超出部分就保持原样,不降噪就好了
// int sub = this->wav_size - (frame_num - 1) * frame_shift;
// this->file_size -= sub;
// this->wav_size -= sub;
// this->data_size -= sub;
}
void Wav_File::FFT(float_t param_real[], float_t param_imag[], int16_t param_n) {
float_t treal, timag, ureal, uimag, arg;
int32_t m, k, j, t, index1, index2;
auto *wreal = (float_t *) malloc(frame_shift * sizeof(float_t));
auto *wimag = (float_t *) malloc(frame_shift * sizeof(float_t));
this->bitrp(param_real, param_imag, param_n);
arg = -2 * PI / param_n;
treal = cosf(arg);
timag = sinf(arg);
wreal[0] = 1.0;
wimag[0] = 0.0;
for (j = 1; j < param_n / 2; j++) {
wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
}
for (m = 2; m <= param_n; m *= 2) {
for (k = 0; k < param_n; k += m) {
for (j = 0; j < m / 2; j++) {
index1 = k + j;
index2 = index1 + m / 2;
t = param_n * j / m;
treal = wreal[t] * param_real[index2] - wimag[t] * param_imag[index2];
timag = wreal[t] * param_imag[index2] + wimag[t] * param_real[index2];
ureal = param_real[index1];
uimag = param_imag[index1];
param_real[index1] = ureal + treal;
param_imag[index1] = uimag + timag;
param_real[index2] = ureal - treal;
param_imag[index2] = uimag - timag;
}
}
}
free(wreal);
free(wimag);
}
void Wav_File::IFFT(float_t param_real[], float_t param_imag[], int16_t param_n) {
float_t treal, timag, ureal, uimag, arg;
int m, k, j, t, index1, index2;
auto *wreal = (float_t *) malloc(frame_shift * sizeof(float_t));
auto *wimag = (float_t *) malloc(frame_shift * sizeof(float_t));
bitrp(param_real, imag, param_n);
arg = 2 * PI / param_n;
treal = cosf(arg);
timag = sinf(arg);
wreal[0] = 1.0;
wimag[0] = 0.0;
for (j = 1; j < param_n / 2; j++) {
wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
}
for (m = 2; m <= param_n; m *= 2) {
for (k = 0; k < param_n; k += m) {
for (j = 0; j < m / 2; j++) {
index1 = k + j;
index2 = index1 + m / 2;
t = param_n * j / m;
treal = wreal[t] * param_real[index2] - wimag[t] * param_imag[index2];
timag = wreal[t] * param_imag[index2] + wimag[t] * param_real[index2];
ureal = param_real[index1];
uimag = param_imag[index1];
param_real[index1] = ureal + treal;
param_imag[index1] = uimag + timag;
param_real[index2] = ureal - treal;
param_imag[index2] = uimag - timag;
}
}
}
for (j = 0; j < param_n; j++) {
real[j] /= param_n;
imag[j] /= param_n;
}
free(wreal);
free(wimag);
}
void Wav_File::bitrp(float_t param_real[], float_t param_imag[], int16_t param_n) {
int i, j, a, b, p;
for (i = 1, p = 0; i < param_n; i *= 2) {
p++;
}
for (i = 0; i < param_n; i++) {
a = i;
b = 0;
for (j = 0; j < p; j++) {
b = (b << 1) + (a & 1); // b = b * 2 + a % 2;
a >>= 1; // a = a / 2;
}
if (b > i) {
std::swap(param_real[i], param_real[b]);
std::swap(param_imag[i], param_imag[b]);
}
}
}
2.3 nfft
不放了 看github吧
3. 其他相关操作
3.1 能量归一化
float_t Wav_File::getNorm(float_t scale) {
if (scale == 0) {
float_t square = 0.0;
for (int i = 0; i < this->wav_size; i++) {
square += powf(this->data[i], 2);
}
scale = sqrtf((float_t) this->wav_size / square);
for (int i = 0; i < this->wav_size; i++) {
this->data[i] *= scale;
}
} else {
for (int i = 0; i < this->wav_size; i++) {
this->data[i] /= scale;
}
}
return scale;
}
3.2 求幅度谱
void Wav_File::getMagnitude() {
for (int i = 0; i < this->frame_num; i++) {
for (int j = 0; j < this->fft_size / 2 + 1; j++) {
this->spec_mag[i][j] = sqrtf(powf(spec_real[i][j], 2) + powf(spec_imag[i][j], 2));
}
}
}
3.3 求相位谱
void Wav_File::getPhase() {
for (int i = 0; i < this->frame_num; i++) {
for (int j = 0; j < this->fft_size / 2 + 1; j++) {
this->spec_pha[i][j] = atan2f(spec_imag[i][j], spec_real[i][j]);
}
}
}
3.4 幅度谱相位谱还原频谱
void Wav_File::magToSpec() {
for (int i = 0; i < this->frame_num; i++) {
for (int j = 0; j < this->fft_size / 2 + 1; j++) {
spec_real[i][j] = this->spec_mag[i][j] * cosf(this->spec_pha[i][j]);
spec_imag[i][j] = this->spec_mag[i][j] * sinf(this->spec_pha[i][j]);
}
}
}
4. 参考链接
[1] WAV文件格式详解
[2] GITHUB源码