Mercurial > hg > lowfreq

changeset 5:bf742ae09443

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Frequency range, and resample to higher frequency to avoid transition band -- and some comments on what is still wrong with it
author Chris Cannam
date Mon, 10 Mar 2014 16:03:53 +0000
parents 9867e53a2592
children e502cf649389
files LowFreq.cpp LowFreq.h
diffstat 2 files changed, 126 insertions(+), 49 deletions(-) [+]
line wrap: on
line diff
--- a/LowFreq.cpp	Fri Mar 07 14:17:47 2014 +0000
+++ b/LowFreq.cpp	Mon Mar 10 16:03:53 2014 +0000
@@ -11,9 +11,13 @@
 using std::endl;
 using std::vector;
 
-static float minP = 0.01;
-static float maxP = 10;
-static float defaultP = 0.1;
+static float minFmin = 0;
+static float maxFmin = 500;
+static float defaultFmin = 0.1;
+
+static float minFmax = 0.1;
+static float maxFmax = 500;
+static float defaultFmax = 10;
 
 static int minN = 1;
 static int maxN = 16384;
@@ -26,7 +30,8 @@
 
 LowFreq::LowFreq(float inputSampleRate) :
     Plugin(inputSampleRate),
-    m_p(defaultP),
+    m_fmin(defaultFmin),
+    m_fmax(defaultFmax),
     m_n(defaultN),
     m_overlap(defaultOverlap),
     m_blockSize(0),
@@ -116,19 +121,29 @@
     ParameterList list;
 
     ParameterDescriptor d;
-    d.identifier = "p";
-    d.name = "Shortest period";
-    d.description = "Period in seconds of the highest-frequency component to include in the spectrogram. That is, 1/f where f is the highest frequency (in Hz) spanned by the spectrogram.";
-    d.unit = "s";
-    d.minValue = minP;
-    d.maxValue = maxP;
-    d.defaultValue = defaultP;
+    d.identifier = "fmin";
+    d.name = "Minimum frequency";
+    d.description = "Low bound for frequencies to include in the spectrogram. This will be rounded up to the next highest bin frequency based on the highest frequency and number of bins.";
+    d.unit = "Hz";
+    d.minValue = minFmin;
+    d.maxValue = maxFmin;
+    d.defaultValue = defaultFmin;
+    d.isQuantized = false;
+    list.push_back(d);
+
+    d.identifier = "fmax";
+    d.name = "Maximum frequency";
+    d.description = "Highest frequency component to include in the spectrogram.";
+    d.unit = "Hz";
+    d.minValue = minFmax;
+    d.maxValue = maxFmax;
+    d.defaultValue = defaultFmax;
     d.isQuantized = false;
     list.push_back(d);
 
     d.identifier = "n";
     d.name = "Number of bins";
-    d.description = "Number of spectrogram bins to calculate.";
+    d.description = "Number of spectrogram bins to return.";
     d.unit = "";
     d.minValue = minN;
     d.maxValue = maxN;
@@ -154,8 +169,10 @@
 float
 LowFreq::getParameter(string identifier) const
 {
-    if (identifier == "p") {
-	return m_p;
+    if (identifier == "fmin") {
+	return m_fmin;
+    } else if (identifier == "fmax") {
+	return m_fmax;
     } else if (identifier == "n") {
 	return m_n;
     } else if (identifier == "overlap") {
@@ -167,8 +184,10 @@
 void
 LowFreq::setParameter(string identifier, float value) 
 {
-    if (identifier == "p") {
-	m_p = value;
+    if (identifier == "fmin") {
+	m_fmin = value;
+    } else if (identifier == "fmax") {
+	m_fmax = value;
     } else if (identifier == "n") {
 	m_n = int(value + 0.01);
     } else if (identifier == "overlap") {
@@ -194,18 +213,6 @@
 {
 }
 
-int
-LowFreq::getTargetSampleRate() const
-{
-    return int(ceil(2.0 / m_p));
-}
-
-int
-LowFreq::getTargetStepSize() const
-{
-    return int(round(m_n * (1.0 - (m_overlap / 100.0))));
-}    
-
 LowFreq::OutputList
 LowFreq::getOutputDescriptors() const
 {
@@ -227,7 +234,7 @@
     
     char namebuf[50];
     for (int i = 0; i <= m_n/2; ++i) {
-	sprintf(namebuf, "%.3gHz", double(getTargetSampleRate()) / m_n * i);
+	sprintf(namebuf, "%.3gHz", double(getOutputBinFrequency(i)));
 	d.binNames.push_back(namebuf);
     }
 
@@ -265,10 +272,17 @@
 	return false;
     }
 
-    if (m_p < minP || m_p > maxP) {
-	cerr << "LowFreq::initialise: ERROR: shortest period " << m_p
-	     << " out of acceptable range " << minP
-	     << " -> " << maxP << endl;
+    if (m_fmin < minFmin || m_fmin > maxFmin) {
+	cerr << "LowFreq::initialise: ERROR: lowest frequency " << m_fmin
+	     << " out of acceptable range " << minFmin
+	     << " -> " << maxFmin << endl;
+	return false;
+    }
+
+    if (m_fmax < minFmax || m_fmax > maxFmax) {
+	cerr << "LowFreq::initialise: ERROR: highest frequency " << m_fmax
+	     << " out of acceptable range " << minFmax
+	     << " -> " << maxFmax << endl;
 	return false;
     }
 
@@ -299,24 +313,79 @@
     delete m_fft;
     delete m_window;
 
-    // We have shortest period p and number of bins n.  The bin
-    // frequency for the highest bin of a spectrogram is fs/2 where fs
-    // is the sample rate of the signal. So we must downsample to a
-    // sample rate of (1/p)*2 and then run an n-point FFT.
+    // We want to return a number of bins n between frequencies fmin
+    // and fmax.
+
+    // The bin frequency for bin n of an n-point fft is fs (the sample
+    // rate), with the top half aliasing the bottom half. So in theory
+    // we can downsample to a sample rate of fmax*2 and then run an
+    // m-point FFT, where m is sufficient to give us n bins remaining
+    // between fmin and fmax. For m then we calculate fmax /
+    // ((fmax-fmin)/n).
+
+    // However, our resampler introduces artifacts close to its cutoff
+    // frequency within the filter transition band so we probably
+    // shouldn't aim to resample to exactly fmax*2. Instead we
+    // resample to fmax*4 and carry out an FFT of twice the expected
+    // number of bins.
 
     m_resampler = new Resampler(int(round(m_inputSampleRate)),
 				getTargetSampleRate());
-    m_fft = new FFT(m_n);
-    m_window = new Window<double>(HanningWindow, m_n);
+    m_fft = new FFT(getFFTSize());
+    m_window = new Window<double>(HanningWindow, getFFTSize());
     m_buffer = std::vector<double>();
 
-    cerr << "LowFreq::reset: min period " << m_p << " so max freq " << 1.0/m_p
-	 << endl;
+    cerr << "LowFreq::reset: freq range " << m_fmin << " to " << m_fmax << endl;
     cerr << "LowFreq::reset: block size " << m_blockSize
 	 << ", input sample rate " << m_inputSampleRate << endl;
     cerr << "LowFreq::reset: target rate " << getTargetSampleRate()
+	 << ", target step " << getTargetStepSize()
 	 << ", resampler latency " << m_resampler->getLatency()
-	 << ", fft size " << m_n << endl;
+	 << ", fft size " << getFFTSize() << endl;
+
+    //!!! output has wrong number of bins
+
+    //!!! not handling resampler latency
+
+    //!!! how do we know the fft size is long enough to capture the low frequencies we want?
+}
+
+int
+LowFreq::getTargetSampleRate() const
+{
+    // See notes in reset()
+    return int(ceil(m_fmax)) * 4;
+}
+
+int
+LowFreq::getTargetStepSize() const
+{
+    return int(round(m_n * (1.0 - (m_overlap / 100.0))));
+}    
+
+int
+LowFreq::getFFTSize() const
+{
+    int fftSize = 2 * int(ceil(m_fmax / ((m_fmax - m_fmin) / m_n)));
+    cerr << "LowFreq::getFFTSize: for range " << m_fmin << " -> " << m_fmax
+	 << " and n = " << m_n << ", fft size is " << fftSize << endl;
+    return fftSize;
+}
+
+int
+LowFreq::getFirstOutputBin() const
+{
+    int first = int(floor(m_fmin / getFFTSize()));
+    cerr << "LowFreq::getFirstOutputBin: for range " << m_fmin << " -> " << m_fmax
+	 << " and n = " << m_n << ", bin is " << first << endl;
+    return first;
+}
+
+float
+LowFreq::getOutputBinFrequency(int i) const
+{
+    return (float(getTargetSampleRate()) / getFFTSize())
+	* (i + getFirstOutputBin());
 }
 
 LowFreq::FeatureSet
@@ -334,7 +403,7 @@
 
     FeatureSet fs;
 
-    while (int(m_buffer.size()) >= m_n) {
+    while (int(m_buffer.size()) >= getFFTSize()) {
 	Feature f = processColumn();
 	fs[0].push_back(f);
 	advance();
@@ -352,7 +421,7 @@
 	m_buffer.push_back(0.0);
     }
 
-    while (int(m_buffer.size()) >= m_n) {
+    while (int(m_buffer.size()) >= getFFTSize()) {
 	Feature f = processColumn();
 	fs[0].push_back(f);
 	advance();
@@ -365,17 +434,21 @@
 LowFreq::processColumn()
 {
     Feature f;
+    
+    int sz = getFFTSize();
 
-    double *realOut = new double[m_n];
-    double *imagOut = new double[m_n];
+    double *realOut = new double[sz];
+    double *imagOut = new double[sz];
 
-    double *windowed = new double[m_n];
+    double *windowed = new double[sz];
     m_window->cut(m_buffer.data(), windowed);	
 
     m_fft->process(false, windowed, 0, realOut, imagOut);
 
     for (int i = 0; i <= m_n/2; ++i) {
-	f.values.push_back(realOut[i] * realOut[i] + imagOut[i] * imagOut[i]);
+	int ix = getFirstOutputBin() + i;
+	f.values.push_back(realOut[ix] * realOut[ix] +
+			   imagOut[ix] * imagOut[ix]);
     }
 
     return f;
--- a/LowFreq.h	Fri Mar 07 14:17:47 2014 +0000
+++ b/LowFreq.h	Mon Mar 10 16:03:53 2014 +0000
@@ -53,8 +53,12 @@
 
     int getTargetSampleRate() const;
     int getTargetStepSize() const;
+    int getFFTSize() const;
+    int getFirstOutputBin() const;
+    float getOutputBinFrequency(int i) const;
 
-    float m_p;
+    float m_fmin;
+    float m_fmax;
     int m_n;
     float m_overlap;