i'm trying use c++ recreate spectrogram function used matlab. function uses short time fourier transform (stft). found c++ code here performs stft. code seems work frequencies want few. found this post similar question following answer:
just take inner product of data complex exponential @ frequency of interest. if g data, substitute f value of frequency want (e.g., 1, 3, 10, ...)
having no background in mathematics, can't figure out how this. inner product part seems simple enough wikipedia page have absolutely no idea means (with regard formula dft)
a complex exponential @ frequency of interest
could explain how might able this? data structure after stft matrix filled complex numbers. don't know how extract desired frequencies.
relevant function, window hamming, , vector of desired frequencies isn't yet input because don't know them:
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 stop reading in chunks? bool shouldstop = false; int numchunkscompleted = 0; int i; // process each chunk of signal while (chunkposition < signallength && !shouldstop) { // copy chunk our buffer (i = 0; < windowsize; i++) { readindex = chunkposition + i; if (readindex < signallength) { // note windowing! data[i][0] = signal[readindex] * window[i]; data[i][1] = 0.0; } else { // have read beyond signal, zero-pad it! data[i][0] = 0.0; data[i][1] = 0.0; shouldstop = true; } } // perform fft on our chunk fftw_execute(plan_forward); // copy first (windowsize/2 + 1) data points spectrogram. // because fft output mirrored nyquist // frequency, second half of data redundant. how // matlab's spectrogram routine works. (i = 0; < 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 formatting, while ends here. return results; }
look goertzel algorithm or filter example code uses computational equivalent of inner product against complex exponential measure presence or magnitude of specific stationary sinusoidal frequency in signal. performance or resolution depend on length of filter , signal.
Comments
Post a Comment