在短时傅里叶变换中获取特定频率的值
Getting values for specific frequencies in a short time fourier transform
我正在尝试使用 C++ 为具有以下答案的类似问题重新创建 spectrogram function used by Matlab. The function uses a Short Time Fourier Transform (STFT). I found some C++ code here that performs a STFT. The code seems to work perfectly for all frequencies but I only want a few. I found this post:
Just take the inner product of your data with a complex exponential at
the frequency of interest. If g is your data, then just substitute for
f the value of the frequency you want (e.g., 1, 3, 10, ...)
没有数学背景,我不知道该怎么做。从 Wikipedia page 看,内积部分看起来很简单,但我完全不知道他的意思(关于 DFT 的公式)
a complex exponential at frequency of interest
有人可以解释一下我是如何做到这一点的吗?我在 STFT 之后的数据结构是一个充满复数的矩阵。我只是不知道如何提取我想要的频率。
相关函数,其中 window
是汉明函数,所需频率的向量还不是输入,因为我不知道如何处理它们:
Matrix<complex<double>> ShortTimeFourierTransform::Calculate(const vector<double> &signal,
const vector<double> &window, int windowSize, int hopSize)
{
int signalLength = signal.size();
int nOverlap = hopSize;
int cols = (signal.size() - nOverlap) / (windowSize - nOverlap);
Matrix<complex<double>> results(window.size(), cols);
int chunkPosition = 0;
int readIndex;
// Should we stop reading in chunks?
bool shouldStop = false;
int numChunksCompleted = 0;
int i;
// Process each chunk of the signal
while (chunkPosition < signalLength && !shouldStop)
{
// Copy the chunk into our buffer
for (i = 0; i < windowSize; i++)
{
readIndex = chunkPosition + i;
if (readIndex < signalLength)
{
// Note the windowing!
data[i][0] = signal[readIndex] * window[i];
data[i][1] = 0.0;
}
else
{
// we have read beyond the signal, so zero-pad it!
data[i][0] = 0.0;
data[i][1] = 0.0;
shouldStop = true;
}
}
// Perform the FFT on our chunk
fftw_execute(plan_forward);
// Copy the first (windowSize/2 + 1) data points into your spectrogram.
// We do this because the FFT output is mirrored about the nyquist
// frequency, so the second half of the data is redundant. This is how
// Matlab's spectrogram routine works.
for (i = 0; i < windowSize / 2 + 1; i++)
{
double real = fft_result[i][0];
double imaginary = fft_result[i][1];
results(i, numChunksCompleted) = complex<double>(real, imaginary);
}
chunkPosition += hopSize;
numChunksCompleted++;
} // Excuse the formatting, the while ends here.
return results;
}
查找 Goertzel algorithm or filter 示例代码,该代码使用内积与复指数的计算等效项来测量信号中特定固定正弦频率的存在或幅度。性能或分辨率将取决于滤波器的长度和您的信号。
我正在尝试使用 C++ 为具有以下答案的类似问题重新创建 spectrogram function used by Matlab. The function uses a Short Time Fourier Transform (STFT). I found some C++ code here that performs a STFT. The code seems to work perfectly for all frequencies but I only want a few. I found this post:
Just take the inner product of your data with a complex exponential at the frequency of interest. If g is your data, then just substitute for f the value of the frequency you want (e.g., 1, 3, 10, ...)
没有数学背景,我不知道该怎么做。从 Wikipedia page 看,内积部分看起来很简单,但我完全不知道他的意思(关于 DFT 的公式)
a complex exponential at frequency of interest
有人可以解释一下我是如何做到这一点的吗?我在 STFT 之后的数据结构是一个充满复数的矩阵。我只是不知道如何提取我想要的频率。
相关函数,其中 window
是汉明函数,所需频率的向量还不是输入,因为我不知道如何处理它们:
Matrix<complex<double>> ShortTimeFourierTransform::Calculate(const vector<double> &signal,
const vector<double> &window, int windowSize, int hopSize)
{
int signalLength = signal.size();
int nOverlap = hopSize;
int cols = (signal.size() - nOverlap) / (windowSize - nOverlap);
Matrix<complex<double>> results(window.size(), cols);
int chunkPosition = 0;
int readIndex;
// Should we stop reading in chunks?
bool shouldStop = false;
int numChunksCompleted = 0;
int i;
// Process each chunk of the signal
while (chunkPosition < signalLength && !shouldStop)
{
// Copy the chunk into our buffer
for (i = 0; i < windowSize; i++)
{
readIndex = chunkPosition + i;
if (readIndex < signalLength)
{
// Note the windowing!
data[i][0] = signal[readIndex] * window[i];
data[i][1] = 0.0;
}
else
{
// we have read beyond the signal, so zero-pad it!
data[i][0] = 0.0;
data[i][1] = 0.0;
shouldStop = true;
}
}
// Perform the FFT on our chunk
fftw_execute(plan_forward);
// Copy the first (windowSize/2 + 1) data points into your spectrogram.
// We do this because the FFT output is mirrored about the nyquist
// frequency, so the second half of the data is redundant. This is how
// Matlab's spectrogram routine works.
for (i = 0; i < windowSize / 2 + 1; i++)
{
double real = fft_result[i][0];
double imaginary = fft_result[i][1];
results(i, numChunksCompleted) = complex<double>(real, imaginary);
}
chunkPosition += hopSize;
numChunksCompleted++;
} // Excuse the formatting, the while ends here.
return results;
}
查找 Goertzel algorithm or filter 示例代码,该代码使用内积与复指数的计算等效项来测量信号中特定固定正弦频率的存在或幅度。性能或分辨率将取决于滤波器的长度和您的信号。