使用 FFT 的移频器
Frequency Shifter Using FFT
我一直在使用 Rosetta Code 提供的原始 FFT 算法开发移频器。我知道要对样本信号进行频移,可以对原始音频应用 FFT,将每个所得正弦波的频率乘以频移因子(用户定义),然后将正弦波加回一起。当我 运行 我的算法时,输出质量极低,就好像算法中没有收集到足够的正弦波来首先接近正确地再现信号一样。该算法在头文件中的 class 中实现,并在其他地方(正确地)调用。
#include <complex>
#include <valarray>
typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;
class FrequencyShifter {
float sampleRate;
public:
FrequencyShifter() {
}
void setSampleRate(float inSampleRate) {
sampleRate = inSampleRate;
}
double abs(double in0) {
if (in0>=0) return in0;
else return -in0;
}
void fft(CArray& x)
{
const size_t N = x.size();
if (N <= 1) return;
// divide
CArray even = x[std::slice(0, N/2, 2)];
CArray odd = x[std::slice(1, N/2, 2)];
// conquer
fft(even);
fft(odd);
// combine
for (size_t k = 0; k < N/2; ++k)
{
Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
x[k ] = even[k] + t;
x[k+N/2] = even[k] - t;
}
}
double convertToReal(double im, double re) {
return sqrt(abs(im*im - re*re));
}
void processBlock(float *inBlock, const int inFramesToProcess, float scale) {
//inFramesToProcess is the amount of samples in inBlock
Complex *copy = new Complex[inFramesToProcess];
for (int frame = 0; frame<inFramesToProcess; frame++) {
copy[frame] = Complex((double)inBlock[frame], 0.0);
}
CArray data(copy, inFramesToProcess);
fft(data);
const float freqoffsets = sampleRate/inFramesToProcess;
for (float x = 0; x<data.size()/2; x++) {
for (float frame = 0; frame<inFramesToProcess; frame++) {
inBlock[(int)frame] = (float)(convertToReal(data[(int)x].imag(), data[(int)x].real())*sin(freqoffsets*x*frame*scale));
}
}
}
};
我假设部分问题是我只包括 sampleRate/inFramesToProcess
正弦波覆盖的频率。发送更大的音频文件(因此更大的 *inBlock
s 和 inFramesToProcess
s)是否会使音频不那么粗糙? 如果不 仅更改参数的值或长度,我将如何完成此操作?
这是 processBlock
的更新版本,其中包含实现频移所需的一些调整,我将在下面进行描述:
void processBlock(float *inBlock, const int inFramesToProcess, float scale) {
//inFramesToProcess is the amount of samples in inBlock
Complex *copy = new Complex[inFramesToProcess];
for (int frame = 0; frame<inFramesToProcess; frame++) {
copy[frame] = Complex((double)inBlock[frame], 0.0);
}
CArray data(copy, inFramesToProcess);
fft(data);
const float freqoffsets = 2.0*PI/inFramesToProcess;
const float normfactor = 2.0/inFramesToProcess;
for (int frame = 0; frame<inFramesToProcess; frame++) {
inBlock[frame] = 0.5*data[0].real();
for (int x = 1; x<data.size()/2; x++) {
float arg = freqoffsets*x*frame*scale;
inBlock[frame] += data[x].real()*cos(arg) - data[x].imag()*sin(arg);
}
inBlock[frame] *= normfactor;
}
}
推导
您从 FFT 获得的频谱是 complex-valued,可以将其视为以正弦波和余弦波表示您的信号。可以使用逆变换来重建 time-domain 波形,逆变换由以下关系给出:
利用频谱对称性,可以表示为:
或等同于:
您可能已经注意到索引 0
和 N/2
处的项是在频域中具有纯实系数的特例。为简单起见,假设频谱没有一直到 N/2
,您可以删除 N/2
项并仍然得到合理的近似值。对于其他条款,您将获得可以实现为
的贡献
normfactor = 2.0/inFramesToProcess;
normfactor*(data[x].real()*cos(arg) - data[x].imag()*sin(arg))
您当然需要将所有这些贡献添加到最终缓冲区 inBlock[frame]
,而不是简单地覆盖以前的结果:
inBlock[frame] += normfactor*(data[x].real()*cos(arg) - data[x].imag()*sin(arg));
// ^^
注意,可以对循环后的最终结果进行归一化,以减少乘法次数。这样做时,我们必须特别注意索引 0 处的 DC 项(其系数为 1/N
而不是 2/N
):
inBlock[frame] = 0.5*data[0].real();
for (int x = 1; x<data.size()/2; x++) {
float arg = freqoffsets*x*frame*scale;
inBlock[frame] += data[x].real()*cos(arg) - data[x].imag()*sin(arg);
}
inBlock[frame] *= normfactor;
最后,在生成音调时,相位参数 arg
到 sin
和 cos
应该是 2*pi*k*n/inFramesToProcess
的形式(在应用 scale
因子),其中 n
是 time-domain 样本索引,k
是频域索引。最终结果是计算出的频率增量 freqoffsets
实际上应该是 2.0*PI/inFramesToProcess
。
备注
- FFT 算法的工作假设是您的基础 time-domain 信号与您的块长度周期是周期性的。因此,块之间可能会出现不连续的声音。
- 未来的读者应该知道,这不会将频谱移动一个常数,而是在频率按乘法因子缩放时挤压或扩展频谱。例如,包含 100-200Hz 分量的信号可能会被压缩到 75-150Hz 0.75 倍。请注意下限是如何下移 25Hz 而上限是下移 50Hz。
我一直在使用 Rosetta Code 提供的原始 FFT 算法开发移频器。我知道要对样本信号进行频移,可以对原始音频应用 FFT,将每个所得正弦波的频率乘以频移因子(用户定义),然后将正弦波加回一起。当我 运行 我的算法时,输出质量极低,就好像算法中没有收集到足够的正弦波来首先接近正确地再现信号一样。该算法在头文件中的 class 中实现,并在其他地方(正确地)调用。
#include <complex>
#include <valarray>
typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;
class FrequencyShifter {
float sampleRate;
public:
FrequencyShifter() {
}
void setSampleRate(float inSampleRate) {
sampleRate = inSampleRate;
}
double abs(double in0) {
if (in0>=0) return in0;
else return -in0;
}
void fft(CArray& x)
{
const size_t N = x.size();
if (N <= 1) return;
// divide
CArray even = x[std::slice(0, N/2, 2)];
CArray odd = x[std::slice(1, N/2, 2)];
// conquer
fft(even);
fft(odd);
// combine
for (size_t k = 0; k < N/2; ++k)
{
Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
x[k ] = even[k] + t;
x[k+N/2] = even[k] - t;
}
}
double convertToReal(double im, double re) {
return sqrt(abs(im*im - re*re));
}
void processBlock(float *inBlock, const int inFramesToProcess, float scale) {
//inFramesToProcess is the amount of samples in inBlock
Complex *copy = new Complex[inFramesToProcess];
for (int frame = 0; frame<inFramesToProcess; frame++) {
copy[frame] = Complex((double)inBlock[frame], 0.0);
}
CArray data(copy, inFramesToProcess);
fft(data);
const float freqoffsets = sampleRate/inFramesToProcess;
for (float x = 0; x<data.size()/2; x++) {
for (float frame = 0; frame<inFramesToProcess; frame++) {
inBlock[(int)frame] = (float)(convertToReal(data[(int)x].imag(), data[(int)x].real())*sin(freqoffsets*x*frame*scale));
}
}
}
};
我假设部分问题是我只包括 sampleRate/inFramesToProcess
正弦波覆盖的频率。发送更大的音频文件(因此更大的 *inBlock
s 和 inFramesToProcess
s)是否会使音频不那么粗糙? 如果不 仅更改参数的值或长度,我将如何完成此操作?
这是 processBlock
的更新版本,其中包含实现频移所需的一些调整,我将在下面进行描述:
void processBlock(float *inBlock, const int inFramesToProcess, float scale) {
//inFramesToProcess is the amount of samples in inBlock
Complex *copy = new Complex[inFramesToProcess];
for (int frame = 0; frame<inFramesToProcess; frame++) {
copy[frame] = Complex((double)inBlock[frame], 0.0);
}
CArray data(copy, inFramesToProcess);
fft(data);
const float freqoffsets = 2.0*PI/inFramesToProcess;
const float normfactor = 2.0/inFramesToProcess;
for (int frame = 0; frame<inFramesToProcess; frame++) {
inBlock[frame] = 0.5*data[0].real();
for (int x = 1; x<data.size()/2; x++) {
float arg = freqoffsets*x*frame*scale;
inBlock[frame] += data[x].real()*cos(arg) - data[x].imag()*sin(arg);
}
inBlock[frame] *= normfactor;
}
}
推导
您从 FFT 获得的频谱是 complex-valued,可以将其视为以正弦波和余弦波表示您的信号。可以使用逆变换来重建 time-domain 波形,逆变换由以下关系给出:
利用频谱对称性,可以表示为:
或等同于:
您可能已经注意到索引 0
和 N/2
处的项是在频域中具有纯实系数的特例。为简单起见,假设频谱没有一直到 N/2
,您可以删除 N/2
项并仍然得到合理的近似值。对于其他条款,您将获得可以实现为
normfactor = 2.0/inFramesToProcess;
normfactor*(data[x].real()*cos(arg) - data[x].imag()*sin(arg))
您当然需要将所有这些贡献添加到最终缓冲区 inBlock[frame]
,而不是简单地覆盖以前的结果:
inBlock[frame] += normfactor*(data[x].real()*cos(arg) - data[x].imag()*sin(arg));
// ^^
注意,可以对循环后的最终结果进行归一化,以减少乘法次数。这样做时,我们必须特别注意索引 0 处的 DC 项(其系数为 1/N
而不是 2/N
):
inBlock[frame] = 0.5*data[0].real();
for (int x = 1; x<data.size()/2; x++) {
float arg = freqoffsets*x*frame*scale;
inBlock[frame] += data[x].real()*cos(arg) - data[x].imag()*sin(arg);
}
inBlock[frame] *= normfactor;
最后,在生成音调时,相位参数 arg
到 sin
和 cos
应该是 2*pi*k*n/inFramesToProcess
的形式(在应用 scale
因子),其中 n
是 time-domain 样本索引,k
是频域索引。最终结果是计算出的频率增量 freqoffsets
实际上应该是 2.0*PI/inFramesToProcess
。
备注
- FFT 算法的工作假设是您的基础 time-domain 信号与您的块长度周期是周期性的。因此,块之间可能会出现不连续的声音。
- 未来的读者应该知道,这不会将频谱移动一个常数,而是在频率按乘法因子缩放时挤压或扩展频谱。例如,包含 100-200Hz 分量的信号可能会被压缩到 75-150Hz 0.75 倍。请注意下限是如何下移 25Hz 而上限是下移 50Hz。