Move SpectrumAnalyst into its own files...
... It's the pure calculation common to the Plot Spectrum window and to spectral editing This removes some dependencies on FreqWindow
This commit is contained in:
Paul Licameli
parent
72f20d9129
commit
6eb0f3aca1
|
|
@ -230,6 +230,8 @@ src/SoundActivatedRecord.cpp
|
|||
src/SoundActivatedRecord.h
|
||||
src/Spectrum.cpp
|
||||
src/Spectrum.h
|
||||
src/SpectrumAnalyst.cpp
|
||||
src/SpectrumAnalyst.h
|
||||
src/SplashDialog.cpp
|
||||
src/SplashDialog.h
|
||||
src/SseMathFuncs.cpp
|
||||
|
|
|
|||
|
|
@ -1313,6 +1313,7 @@
|
|||
5EA0182A1EC7B226001F2996 /* NoteTrackVRulerControls.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 5EA0181F1EC7B226001F2996 /* NoteTrackVRulerControls.cpp */; };
|
||||
5EA0182B1EC7B226001F2996 /* WaveTrackControls.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 5EA018231EC7B226001F2996 /* WaveTrackControls.cpp */; };
|
||||
5EA0182D1EC7B226001F2996 /* WaveTrackVRulerControls.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 5EA018261EC7B226001F2996 /* WaveTrackVRulerControls.cpp */; };
|
||||
5EAF751923C0EA4F00E94479 /* SpectrumAnalyst.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 5EAF751723C0EA4E00E94479 /* SpectrumAnalyst.cpp */; };
|
||||
5EB15A2022A94043009FEC89 /* ProjectHistory.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 5EB15A1E22A94043009FEC89 /* ProjectHistory.cpp */; };
|
||||
5EBD1C9422D11DAF00299FD4 /* SpectrumVZoomHandle.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 5EBD1C9022D11DAF00299FD4 /* SpectrumVZoomHandle.cpp */; };
|
||||
5EBD1C9522D11DAF00299FD4 /* WaveformVZoomHandle.cpp in Sources */ = {isa = PBXBuildFile; fileRef = 5EBD1C9222D11DAF00299FD4 /* WaveformVZoomHandle.cpp */; };
|
||||
|
|
@ -3392,6 +3393,8 @@
|
|||
5EA018241EC7B226001F2996 /* WaveTrackControls.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = WaveTrackControls.h; sourceTree = "<group>"; };
|
||||
5EA018261EC7B226001F2996 /* WaveTrackVRulerControls.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = WaveTrackVRulerControls.cpp; sourceTree = "<group>"; };
|
||||
5EA018271EC7B226001F2996 /* WaveTrackVRulerControls.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = WaveTrackVRulerControls.h; sourceTree = "<group>"; };
|
||||
5EAF751723C0EA4E00E94479 /* SpectrumAnalyst.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = SpectrumAnalyst.cpp; sourceTree = "<group>"; };
|
||||
5EAF751823C0EA4E00E94479 /* SpectrumAnalyst.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = SpectrumAnalyst.h; sourceTree = "<group>"; };
|
||||
5EB15A1E22A94043009FEC89 /* ProjectHistory.cpp */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.cpp.cpp; path = ProjectHistory.cpp; sourceTree = "<group>"; };
|
||||
5EB15A1F22A94043009FEC89 /* ProjectHistory.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = ProjectHistory.h; sourceTree = "<group>"; };
|
||||
5EB9EA281D5B81270050AF40 /* ImportForwards.h */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.c.h; path = ImportForwards.h; sourceTree = "<group>"; };
|
||||
|
|
@ -4644,6 +4647,8 @@
|
|||
2860BA210E0F0D8600A13878 /* SoundActivatedRecord.h */,
|
||||
1790B0DE09883BFD008A330A /* Spectrum.cpp */,
|
||||
1790B0DF09883BFD008A330A /* Spectrum.h */,
|
||||
5EAF751723C0EA4E00E94479 /* SpectrumAnalyst.cpp */,
|
||||
5EAF751823C0EA4E00E94479 /* SpectrumAnalyst.h */,
|
||||
28501E9F0CEECEF80029ABAA /* SplashDialog.cpp */,
|
||||
28501EA00CEECEF80029ABAA /* SplashDialog.h */,
|
||||
EDFCEBA418894B2A00C98E51 /* SseMathFuncs.cpp */,
|
||||
|
|
@ -8713,6 +8718,7 @@
|
|||
1790B16F09883BFD008A330A /* RawAudioGuess.cpp in Sources */,
|
||||
1790B17009883BFD008A330A /* Internat.cpp in Sources */,
|
||||
1790B17109883BFD008A330A /* LabelTrack.cpp in Sources */,
|
||||
5EAF751923C0EA4F00E94479 /* SpectrumAnalyst.cpp in Sources */,
|
||||
1790B17309883BFD008A330A /* LangChoice.cpp in Sources */,
|
||||
5EF3E651203FBAFB006C6882 /* LoadCommands.cpp in Sources */,
|
||||
282B70331B682342009A1618 /* WaveformPrefs.cpp in Sources */,
|
||||
|
|
|
|||
|
|
@ -24,14 +24,6 @@ This class actually does the graph display.
|
|||
Has a feature that finds peaks and reports their value as you move
|
||||
the mouse around.
|
||||
|
||||
*//****************************************************************//**
|
||||
|
||||
\class SpectrumAnalyst
|
||||
\brief Used for finding the peaks, for snapping to peaks.
|
||||
|
||||
This class is used to do the 'find peaks' snapping both in FreqPlot
|
||||
and in the spectrogram spectral selection.
|
||||
|
||||
*//*******************************************************************/
|
||||
|
||||
/*
|
||||
|
|
@ -52,9 +44,9 @@ and in the spectrogram spectral selection.
|
|||
#include <wx/checkbox.h>
|
||||
#include <wx/choice.h>
|
||||
#include <wx/dcclient.h>
|
||||
#include <wx/dcmemory.h>
|
||||
#include <wx/font.h>
|
||||
#include <wx/image.h>
|
||||
#include <wx/dcmemory.h>
|
||||
#include <wx/file.h>
|
||||
#include <wx/filedlg.h>
|
||||
#include <wx/intl.h>
|
||||
|
|
@ -188,17 +180,6 @@ BEGIN_EVENT_TABLE(FrequencyPlotDialog, wxDialogWrapper)
|
|||
EVT_COMMAND(wxID_ANY, EVT_FREQWINDOW_RECALC, FrequencyPlotDialog::OnRecalc)
|
||||
END_EVENT_TABLE()
|
||||
|
||||
SpectrumAnalyst::SpectrumAnalyst()
|
||||
: mAlg(Spectrum)
|
||||
, mRate(0.0)
|
||||
, mWindowSize(0)
|
||||
{
|
||||
}
|
||||
|
||||
SpectrumAnalyst::~SpectrumAnalyst()
|
||||
{
|
||||
}
|
||||
|
||||
FrequencyPlotDialog::FrequencyPlotDialog(wxWindow * parent, wxWindowID id,
|
||||
AudacityProject &project,
|
||||
const TranslatableString & title,
|
||||
|
|
@ -1156,471 +1137,3 @@ void FreqPlot::OnMouseEvent(wxMouseEvent & event)
|
|||
freqWindow->PlotMouseEvent(event);
|
||||
}
|
||||
|
||||
FreqGauge::FreqGauge(wxWindow * parent, wxWindowID winid)
|
||||
: wxStatusBar(parent, winid, wxST_SIZEGRIP)
|
||||
{
|
||||
mRange = 0;
|
||||
}
|
||||
|
||||
void FreqGauge::SetRange(int range, int bar, int gap)
|
||||
{
|
||||
mRange = range;
|
||||
mBar = bar;
|
||||
mGap = gap;
|
||||
|
||||
GetFieldRect(0, mRect);
|
||||
mRect.Inflate(-1);
|
||||
|
||||
mInterval = mRange / (mRect.width / (mBar + mGap));
|
||||
mRect.width = mBar;
|
||||
mMargin = mRect.x;
|
||||
mLast = -1;
|
||||
|
||||
Update();
|
||||
}
|
||||
|
||||
void FreqGauge::SetValue(int value)
|
||||
{
|
||||
mCur = value / mInterval;
|
||||
|
||||
if (mCur != mLast)
|
||||
{
|
||||
wxClientDC dc(this);
|
||||
dc.SetPen(*wxTRANSPARENT_PEN);
|
||||
dc.SetBrush(wxColour(100, 100, 220));
|
||||
|
||||
while (mLast < mCur)
|
||||
{
|
||||
mLast++;
|
||||
mRect.x = mMargin + mLast * (mBar + mGap);
|
||||
dc.DrawRectangle(mRect);
|
||||
}
|
||||
Update();
|
||||
}
|
||||
}
|
||||
|
||||
void FreqGauge::Reset()
|
||||
{
|
||||
mRange = 0;
|
||||
Refresh(true);
|
||||
}
|
||||
|
||||
bool SpectrumAnalyst::Calculate(Algorithm alg, int windowFunc,
|
||||
size_t windowSize, double rate,
|
||||
const float *data, size_t dataLen,
|
||||
float *pYMin, float *pYMax,
|
||||
FreqGauge *progress)
|
||||
{
|
||||
// Wipe old data
|
||||
mProcessed.resize(0);
|
||||
mRate = 0.0;
|
||||
mWindowSize = 0;
|
||||
|
||||
// Validate inputs
|
||||
int f = NumWindowFuncs();
|
||||
|
||||
if (!(windowSize >= 32 && windowSize <= 65536 &&
|
||||
alg >= SpectrumAnalyst::Spectrum &&
|
||||
alg < SpectrumAnalyst::NumAlgorithms &&
|
||||
windowFunc >= 0 && windowFunc < f)) {
|
||||
return false;
|
||||
}
|
||||
|
||||
if (dataLen < windowSize) {
|
||||
return false;
|
||||
}
|
||||
|
||||
// Now repopulate
|
||||
mRate = rate;
|
||||
mWindowSize = windowSize;
|
||||
mAlg = alg;
|
||||
|
||||
auto half = mWindowSize / 2;
|
||||
mProcessed.resize(mWindowSize);
|
||||
|
||||
Floats in{ mWindowSize };
|
||||
Floats out{ mWindowSize };
|
||||
Floats out2{ mWindowSize };
|
||||
Floats win{ mWindowSize };
|
||||
|
||||
for (size_t i = 0; i < mWindowSize; i++) {
|
||||
mProcessed[i] = 0.0f;
|
||||
win[i] = 1.0f;
|
||||
}
|
||||
|
||||
WindowFunc(windowFunc, mWindowSize, win.get());
|
||||
|
||||
// Scale window such that an amplitude of 1.0 in the time domain
|
||||
// shows an amplitude of 0dB in the frequency domain
|
||||
double wss = 0;
|
||||
for (size_t i = 0; i<mWindowSize; i++)
|
||||
wss += win[i];
|
||||
if(wss > 0)
|
||||
wss = 4.0 / (wss*wss);
|
||||
else
|
||||
wss = 1.0;
|
||||
|
||||
if (progress) {
|
||||
progress->SetRange(dataLen);
|
||||
}
|
||||
|
||||
size_t start = 0;
|
||||
int windows = 0;
|
||||
while (start + mWindowSize <= dataLen) {
|
||||
for (size_t i = 0; i < mWindowSize; i++)
|
||||
in[i] = win[i] * data[start + i];
|
||||
|
||||
switch (alg) {
|
||||
case Spectrum:
|
||||
PowerSpectrum(mWindowSize, in.get(), out.get());
|
||||
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] += out[i];
|
||||
break;
|
||||
|
||||
case Autocorrelation:
|
||||
case CubeRootAutocorrelation:
|
||||
case EnhancedAutocorrelation:
|
||||
|
||||
// Take FFT
|
||||
RealFFT(mWindowSize, in.get(), out.get(), out2.get());
|
||||
// Compute power
|
||||
for (size_t i = 0; i < mWindowSize; i++)
|
||||
in[i] = (out[i] * out[i]) + (out2[i] * out2[i]);
|
||||
|
||||
if (alg == Autocorrelation) {
|
||||
for (size_t i = 0; i < mWindowSize; i++)
|
||||
in[i] = sqrt(in[i]);
|
||||
}
|
||||
if (alg == CubeRootAutocorrelation ||
|
||||
alg == EnhancedAutocorrelation) {
|
||||
// Tolonen and Karjalainen recommend taking the cube root
|
||||
// of the power, instead of the square root
|
||||
|
||||
for (size_t i = 0; i < mWindowSize; i++)
|
||||
in[i] = pow(in[i], 1.0f / 3.0f);
|
||||
}
|
||||
// Take FFT
|
||||
RealFFT(mWindowSize, in.get(), out.get(), out2.get());
|
||||
|
||||
// Take real part of result
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] += out[i];
|
||||
break;
|
||||
|
||||
case Cepstrum:
|
||||
RealFFT(mWindowSize, in.get(), out.get(), out2.get());
|
||||
|
||||
// Compute log power
|
||||
// Set a sane lower limit assuming maximum time amplitude of 1.0
|
||||
{
|
||||
float power;
|
||||
float minpower = 1e-20*mWindowSize*mWindowSize;
|
||||
for (size_t i = 0; i < mWindowSize; i++)
|
||||
{
|
||||
power = (out[i] * out[i]) + (out2[i] * out2[i]);
|
||||
if(power < minpower)
|
||||
in[i] = log(minpower);
|
||||
else
|
||||
in[i] = log(power);
|
||||
}
|
||||
// Take IFFT
|
||||
InverseRealFFT(mWindowSize, in.get(), NULL, out.get());
|
||||
|
||||
// Take real part of result
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] += out[i];
|
||||
}
|
||||
|
||||
break;
|
||||
|
||||
default:
|
||||
wxASSERT(false);
|
||||
break;
|
||||
} //switch
|
||||
|
||||
// Update the progress bar
|
||||
if (progress) {
|
||||
progress->SetValue(start);
|
||||
}
|
||||
|
||||
start += half;
|
||||
windows++;
|
||||
}
|
||||
|
||||
if (progress) {
|
||||
// Reset for next time
|
||||
progress->Reset();
|
||||
}
|
||||
|
||||
float mYMin = 1000000, mYMax = -1000000;
|
||||
double scale;
|
||||
switch (alg) {
|
||||
case Spectrum:
|
||||
// Convert to decibels
|
||||
mYMin = 1000000.;
|
||||
mYMax = -1000000.;
|
||||
scale = wss / (double)windows;
|
||||
for (size_t i = 0; i < half; i++)
|
||||
{
|
||||
mProcessed[i] = 10 * log10(mProcessed[i] * scale);
|
||||
if(mProcessed[i] > mYMax)
|
||||
mYMax = mProcessed[i];
|
||||
else if(mProcessed[i] < mYMin)
|
||||
mYMin = mProcessed[i];
|
||||
}
|
||||
break;
|
||||
|
||||
case Autocorrelation:
|
||||
case CubeRootAutocorrelation:
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] = mProcessed[i] / windows;
|
||||
|
||||
// Find min/max
|
||||
mYMin = mProcessed[0];
|
||||
mYMax = mProcessed[0];
|
||||
for (size_t i = 1; i < half; i++)
|
||||
if (mProcessed[i] > mYMax)
|
||||
mYMax = mProcessed[i];
|
||||
else if (mProcessed[i] < mYMin)
|
||||
mYMin = mProcessed[i];
|
||||
break;
|
||||
|
||||
case EnhancedAutocorrelation:
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] = mProcessed[i] / windows;
|
||||
|
||||
// Peak Pruning as described by Tolonen and Karjalainen, 2000
|
||||
|
||||
// Clip at zero, copy to temp array
|
||||
for (size_t i = 0; i < half; i++) {
|
||||
if (mProcessed[i] < 0.0)
|
||||
mProcessed[i] = float(0.0);
|
||||
out[i] = mProcessed[i];
|
||||
}
|
||||
|
||||
// Subtract a time-doubled signal (linearly interp.) from the original
|
||||
// (clipped) signal
|
||||
for (size_t i = 0; i < half; i++)
|
||||
if ((i % 2) == 0)
|
||||
mProcessed[i] -= out[i / 2];
|
||||
else
|
||||
mProcessed[i] -= ((out[i / 2] + out[i / 2 + 1]) / 2);
|
||||
|
||||
// Clip at zero again
|
||||
for (size_t i = 0; i < half; i++)
|
||||
if (mProcessed[i] < 0.0)
|
||||
mProcessed[i] = float(0.0);
|
||||
|
||||
// Find NEW min/max
|
||||
mYMin = mProcessed[0];
|
||||
mYMax = mProcessed[0];
|
||||
for (size_t i = 1; i < half; i++)
|
||||
if (mProcessed[i] > mYMax)
|
||||
mYMax = mProcessed[i];
|
||||
else if (mProcessed[i] < mYMin)
|
||||
mYMin = mProcessed[i];
|
||||
break;
|
||||
|
||||
case Cepstrum:
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] = mProcessed[i] / windows;
|
||||
|
||||
// Find min/max, ignoring first and last few values
|
||||
{
|
||||
size_t ignore = 4;
|
||||
mYMin = mProcessed[ignore];
|
||||
mYMax = mProcessed[ignore];
|
||||
for (size_t i = ignore + 1; i + ignore < half; i++)
|
||||
if (mProcessed[i] > mYMax)
|
||||
mYMax = mProcessed[i];
|
||||
else if (mProcessed[i] < mYMin)
|
||||
mYMin = mProcessed[i];
|
||||
}
|
||||
break;
|
||||
|
||||
default:
|
||||
wxASSERT(false);
|
||||
break;
|
||||
}
|
||||
|
||||
if (pYMin)
|
||||
*pYMin = mYMin;
|
||||
if (pYMax)
|
||||
*pYMax = mYMax;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
const float *SpectrumAnalyst::GetProcessed() const
|
||||
{
|
||||
return &mProcessed[0];
|
||||
}
|
||||
|
||||
int SpectrumAnalyst::GetProcessedSize() const
|
||||
{
|
||||
return mProcessed.size() / 2;
|
||||
}
|
||||
|
||||
float SpectrumAnalyst::GetProcessedValue(float freq0, float freq1) const
|
||||
{
|
||||
float bin0, bin1, binwidth;
|
||||
|
||||
if (mAlg == Spectrum) {
|
||||
bin0 = freq0 * mWindowSize / mRate;
|
||||
bin1 = freq1 * mWindowSize / mRate;
|
||||
} else {
|
||||
bin0 = freq0 * mRate;
|
||||
bin1 = freq1 * mRate;
|
||||
}
|
||||
binwidth = bin1 - bin0;
|
||||
|
||||
float value = float(0.0);
|
||||
|
||||
if (binwidth < 1.0) {
|
||||
float binmid = (bin0 + bin1) / 2.0;
|
||||
int ibin = (int)(binmid) - 1;
|
||||
if (ibin < 1)
|
||||
ibin = 1;
|
||||
if (ibin >= GetProcessedSize() - 3)
|
||||
ibin = std::max(0, GetProcessedSize() - 4);
|
||||
|
||||
value = CubicInterpolate(mProcessed[ibin],
|
||||
mProcessed[ibin + 1],
|
||||
mProcessed[ibin + 2],
|
||||
mProcessed[ibin + 3], binmid - ibin);
|
||||
|
||||
} else {
|
||||
if (bin0 < 0)
|
||||
bin0 = 0;
|
||||
if (bin1 >= GetProcessedSize())
|
||||
bin1 = GetProcessedSize() - 1;
|
||||
|
||||
if ((int)(bin1) > (int)(bin0))
|
||||
value += mProcessed[(int)(bin0)] * ((int)(bin0) + 1 - bin0);
|
||||
bin0 = 1 + (int)(bin0);
|
||||
while (bin0 < (int)(bin1)) {
|
||||
value += mProcessed[(int)(bin0)];
|
||||
bin0 += 1.0;
|
||||
}
|
||||
value += mProcessed[(int)(bin1)] * (bin1 - (int)(bin1));
|
||||
|
||||
value /= binwidth;
|
||||
}
|
||||
|
||||
return value;
|
||||
}
|
||||
|
||||
float SpectrumAnalyst::FindPeak(float xPos, float *pY) const
|
||||
{
|
||||
float bestpeak = 0.0f;
|
||||
float bestValue = 0.0;
|
||||
if (GetProcessedSize() > 1) {
|
||||
bool up = (mProcessed[1] > mProcessed[0]);
|
||||
float bestdist = 1000000;
|
||||
for (int bin = 3; bin < GetProcessedSize() - 1; bin++) {
|
||||
bool nowUp = mProcessed[bin] > mProcessed[bin - 1];
|
||||
if (!nowUp && up) {
|
||||
// Local maximum. Find actual value by cubic interpolation
|
||||
int leftbin = bin - 2;
|
||||
/*
|
||||
if (leftbin < 1)
|
||||
leftbin = 1;
|
||||
*/
|
||||
float valueAtMax = 0.0;
|
||||
float max = leftbin + CubicMaximize(mProcessed[leftbin],
|
||||
mProcessed[leftbin + 1],
|
||||
mProcessed[leftbin + 2],
|
||||
mProcessed[leftbin + 3],
|
||||
&valueAtMax);
|
||||
|
||||
float thispeak;
|
||||
if (mAlg == Spectrum)
|
||||
thispeak = max * mRate / mWindowSize;
|
||||
else
|
||||
thispeak = max / mRate;
|
||||
|
||||
if (fabs(thispeak - xPos) < bestdist) {
|
||||
bestpeak = thispeak;
|
||||
bestdist = fabs(thispeak - xPos);
|
||||
bestValue = valueAtMax;
|
||||
// Should this test come after the enclosing if?
|
||||
if (thispeak > xPos)
|
||||
break;
|
||||
}
|
||||
}
|
||||
up = nowUp;
|
||||
}
|
||||
}
|
||||
|
||||
if (pY)
|
||||
*pY = bestValue;
|
||||
return bestpeak;
|
||||
}
|
||||
|
||||
// If f(0)=y0, f(1)=y1, f(2)=y2, and f(3)=y3, this function finds
|
||||
// the degree-three polynomial which best fits these points and
|
||||
// returns the value of this polynomial at a value x. Usually
|
||||
// 0 < x < 3
|
||||
|
||||
float SpectrumAnalyst::CubicInterpolate(float y0, float y1, float y2, float y3, float x) const
|
||||
{
|
||||
float a, b, c, d;
|
||||
|
||||
a = y0 / -6.0 + y1 / 2.0 - y2 / 2.0 + y3 / 6.0;
|
||||
b = y0 - 5.0 * y1 / 2.0 + 2.0 * y2 - y3 / 2.0;
|
||||
c = -11.0 * y0 / 6.0 + 3.0 * y1 - 3.0 * y2 / 2.0 + y3 / 3.0;
|
||||
d = y0;
|
||||
|
||||
float xx = x * x;
|
||||
float xxx = xx * x;
|
||||
|
||||
return (a * xxx + b * xx + c * x + d);
|
||||
}
|
||||
|
||||
float SpectrumAnalyst::CubicMaximize(float y0, float y1, float y2, float y3, float * max) const
|
||||
{
|
||||
// Find coefficients of cubic
|
||||
|
||||
float a, b, c, d;
|
||||
|
||||
a = y0 / -6.0 + y1 / 2.0 - y2 / 2.0 + y3 / 6.0;
|
||||
b = y0 - 5.0 * y1 / 2.0 + 2.0 * y2 - y3 / 2.0;
|
||||
c = -11.0 * y0 / 6.0 + 3.0 * y1 - 3.0 * y2 / 2.0 + y3 / 3.0;
|
||||
d = y0;
|
||||
|
||||
// Take derivative
|
||||
|
||||
float da, db, dc;
|
||||
|
||||
da = 3 * a;
|
||||
db = 2 * b;
|
||||
dc = c;
|
||||
|
||||
// Find zeroes of derivative using quadratic equation
|
||||
|
||||
float discriminant = db * db - 4 * da * dc;
|
||||
if (discriminant < 0.0)
|
||||
return float(-1.0); // error
|
||||
|
||||
float x1 = (-db + sqrt(discriminant)) / (2 * da);
|
||||
float x2 = (-db - sqrt(discriminant)) / (2 * da);
|
||||
|
||||
// The one which corresponds to a local _maximum_ in the
|
||||
// cubic is the one we want - the one with a negative
|
||||
// second derivative
|
||||
|
||||
float dda = 2 * da;
|
||||
float ddb = db;
|
||||
|
||||
if (dda * x1 + ddb < 0)
|
||||
{
|
||||
*max = a*x1*x1*x1+b*x1*x1+c*x1+d;
|
||||
return x1;
|
||||
}
|
||||
else
|
||||
{
|
||||
*max = a*x2*x2*x2+b*x2*x2+c*x2+d;
|
||||
return x2;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -15,6 +15,7 @@
|
|||
#include <wx/font.h> // member variable
|
||||
#include <wx/statusbr.h> // to inherit
|
||||
#include "SampleFormat.h"
|
||||
#include "SpectrumAnalyst.h"
|
||||
#include "widgets/wxPanelWrapper.h" // to inherit
|
||||
|
||||
class wxMemoryDC;
|
||||
|
|
@ -32,68 +33,6 @@ class RulerPanel;
|
|||
|
||||
DECLARE_EXPORTED_EVENT_TYPE(AUDACITY_DLL_API, EVT_FREQWINDOW_RECALC, -1);
|
||||
|
||||
class SpectrumAnalyst
|
||||
{
|
||||
public:
|
||||
|
||||
enum Algorithm {
|
||||
Spectrum,
|
||||
Autocorrelation,
|
||||
CubeRootAutocorrelation,
|
||||
EnhancedAutocorrelation,
|
||||
Cepstrum,
|
||||
|
||||
NumAlgorithms
|
||||
};
|
||||
|
||||
SpectrumAnalyst();
|
||||
~SpectrumAnalyst();
|
||||
|
||||
// Return true iff successful
|
||||
bool Calculate(Algorithm alg,
|
||||
int windowFunc, // see FFT.h for values
|
||||
size_t windowSize, double rate,
|
||||
const float *data, size_t dataLen,
|
||||
float *pYMin = NULL, float *pYMax = NULL, // outputs
|
||||
FreqGauge *progress = NULL);
|
||||
|
||||
const float *GetProcessed() const;
|
||||
int GetProcessedSize() const;
|
||||
|
||||
float GetProcessedValue(float freq0, float freq1) const;
|
||||
float FindPeak(float xPos, float *pY) const;
|
||||
|
||||
private:
|
||||
float CubicInterpolate(float y0, float y1, float y2, float y3, float x) const;
|
||||
float CubicMaximize(float y0, float y1, float y2, float y3, float * max) const;
|
||||
|
||||
private:
|
||||
Algorithm mAlg;
|
||||
double mRate;
|
||||
size_t mWindowSize;
|
||||
std::vector<float> mProcessed;
|
||||
};
|
||||
|
||||
class FreqGauge final : public wxStatusBar
|
||||
{
|
||||
public:
|
||||
FreqGauge(wxWindow * parent, wxWindowID winid);
|
||||
|
||||
void SetRange(int range, int bar = 12, int gap = 3);
|
||||
void SetValue(int value);
|
||||
void Reset();
|
||||
|
||||
private:
|
||||
wxRect mRect;
|
||||
int mRange;
|
||||
int mCur;
|
||||
int mLast;
|
||||
int mInterval;
|
||||
int mBar;
|
||||
int mGap;
|
||||
int mMargin;
|
||||
};
|
||||
|
||||
class FreqPlot final : public wxWindow
|
||||
{
|
||||
public:
|
||||
|
|
|
|||
|
|
@ -276,6 +276,8 @@ audacity_SOURCES = \
|
|||
SoundActivatedRecord.h \
|
||||
Spectrum.cpp \
|
||||
Spectrum.h \
|
||||
SpectrumAnalyst.cpp \
|
||||
SpectrumAnalyst.h \
|
||||
SplashDialog.cpp \
|
||||
SplashDialog.h \
|
||||
SseMathFuncs.cpp \
|
||||
|
|
|
|||
|
|
@ -343,7 +343,8 @@ am__audacity_SOURCES_DIST = BlockFile.cpp BlockFile.h DirManager.cpp \
|
|||
ShuttleGetDefinition.cpp ShuttleGetDefinition.h ShuttleGui.cpp \
|
||||
ShuttleGui.h ShuttlePrefs.cpp ShuttlePrefs.h Snap.cpp Snap.h \
|
||||
SoundActivatedRecord.cpp SoundActivatedRecord.h Spectrum.cpp \
|
||||
Spectrum.h SplashDialog.cpp SplashDialog.h SseMathFuncs.cpp \
|
||||
Spectrum.h SpectrumAnalyst.cpp SpectrumAnalyst.h \
|
||||
SplashDialog.cpp SplashDialog.h SseMathFuncs.cpp \
|
||||
SseMathFuncs.h Tags.cpp Tags.h Theme.cpp Theme.h \
|
||||
ThemeAsCeeCode.h TimeDialog.cpp TimeDialog.h \
|
||||
TimerRecordDialog.cpp TimerRecordDialog.h TimeTrack.cpp \
|
||||
|
|
@ -727,7 +728,8 @@ am_audacity_OBJECTS = $(am__objects_1) audacity-AboutDialog.$(OBJEXT) \
|
|||
audacity-ShuttleGui.$(OBJEXT) audacity-ShuttlePrefs.$(OBJEXT) \
|
||||
audacity-Snap.$(OBJEXT) \
|
||||
audacity-SoundActivatedRecord.$(OBJEXT) \
|
||||
audacity-Spectrum.$(OBJEXT) audacity-SplashDialog.$(OBJEXT) \
|
||||
audacity-Spectrum.$(OBJEXT) audacity-SpectrumAnalyst.$(OBJEXT) \
|
||||
audacity-SplashDialog.$(OBJEXT) \
|
||||
audacity-SseMathFuncs.$(OBJEXT) audacity-Tags.$(OBJEXT) \
|
||||
audacity-Theme.$(OBJEXT) audacity-TimeDialog.$(OBJEXT) \
|
||||
audacity-TimerRecordDialog.$(OBJEXT) \
|
||||
|
|
@ -1121,6 +1123,7 @@ am__depfiles_remade = ./$(DEPDIR)/audacity-AColor.Po \
|
|||
./$(DEPDIR)/audacity-Snap.Po \
|
||||
./$(DEPDIR)/audacity-SoundActivatedRecord.Po \
|
||||
./$(DEPDIR)/audacity-Spectrum.Po \
|
||||
./$(DEPDIR)/audacity-SpectrumAnalyst.Po \
|
||||
./$(DEPDIR)/audacity-SplashDialog.Po \
|
||||
./$(DEPDIR)/audacity-SseMathFuncs.Po \
|
||||
./$(DEPDIR)/audacity-Tags.Po ./$(DEPDIR)/audacity-Theme.Po \
|
||||
|
|
@ -1868,7 +1871,8 @@ audacity_SOURCES = $(libaudacity_la_SOURCES) AboutDialog.cpp \
|
|||
ShuttleGetDefinition.cpp ShuttleGetDefinition.h ShuttleGui.cpp \
|
||||
ShuttleGui.h ShuttlePrefs.cpp ShuttlePrefs.h Snap.cpp Snap.h \
|
||||
SoundActivatedRecord.cpp SoundActivatedRecord.h Spectrum.cpp \
|
||||
Spectrum.h SplashDialog.cpp SplashDialog.h SseMathFuncs.cpp \
|
||||
Spectrum.h SpectrumAnalyst.cpp SpectrumAnalyst.h \
|
||||
SplashDialog.cpp SplashDialog.h SseMathFuncs.cpp \
|
||||
SseMathFuncs.h Tags.cpp Tags.h Theme.cpp Theme.h \
|
||||
ThemeAsCeeCode.h TimeDialog.cpp TimeDialog.h \
|
||||
TimerRecordDialog.cpp TimerRecordDialog.h TimeTrack.cpp \
|
||||
|
|
@ -3139,6 +3143,7 @@ distclean-compile:
|
|||
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/audacity-Snap.Po@am__quote@ # am--include-marker
|
||||
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/audacity-SoundActivatedRecord.Po@am__quote@ # am--include-marker
|
||||
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/audacity-Spectrum.Po@am__quote@ # am--include-marker
|
||||
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/audacity-SpectrumAnalyst.Po@am__quote@ # am--include-marker
|
||||
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/audacity-SplashDialog.Po@am__quote@ # am--include-marker
|
||||
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/audacity-SseMathFuncs.Po@am__quote@ # am--include-marker
|
||||
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/audacity-Tags.Po@am__quote@ # am--include-marker
|
||||
|
|
@ -4997,6 +5002,20 @@ audacity-Spectrum.obj: Spectrum.cpp
|
|||
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
|
||||
@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(audacity_CPPFLAGS) $(CPPFLAGS) $(audacity_CXXFLAGS) $(CXXFLAGS) -c -o audacity-Spectrum.obj `if test -f 'Spectrum.cpp'; then $(CYGPATH_W) 'Spectrum.cpp'; else $(CYGPATH_W) '$(srcdir)/Spectrum.cpp'; fi`
|
||||
|
||||
audacity-SpectrumAnalyst.o: SpectrumAnalyst.cpp
|
||||
@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(audacity_CPPFLAGS) $(CPPFLAGS) $(audacity_CXXFLAGS) $(CXXFLAGS) -MT audacity-SpectrumAnalyst.o -MD -MP -MF $(DEPDIR)/audacity-SpectrumAnalyst.Tpo -c -o audacity-SpectrumAnalyst.o `test -f 'SpectrumAnalyst.cpp' || echo '$(srcdir)/'`SpectrumAnalyst.cpp
|
||||
@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/audacity-SpectrumAnalyst.Tpo $(DEPDIR)/audacity-SpectrumAnalyst.Po
|
||||
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='SpectrumAnalyst.cpp' object='audacity-SpectrumAnalyst.o' libtool=no @AMDEPBACKSLASH@
|
||||
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
|
||||
@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(audacity_CPPFLAGS) $(CPPFLAGS) $(audacity_CXXFLAGS) $(CXXFLAGS) -c -o audacity-SpectrumAnalyst.o `test -f 'SpectrumAnalyst.cpp' || echo '$(srcdir)/'`SpectrumAnalyst.cpp
|
||||
|
||||
audacity-SpectrumAnalyst.obj: SpectrumAnalyst.cpp
|
||||
@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(audacity_CPPFLAGS) $(CPPFLAGS) $(audacity_CXXFLAGS) $(CXXFLAGS) -MT audacity-SpectrumAnalyst.obj -MD -MP -MF $(DEPDIR)/audacity-SpectrumAnalyst.Tpo -c -o audacity-SpectrumAnalyst.obj `if test -f 'SpectrumAnalyst.cpp'; then $(CYGPATH_W) 'SpectrumAnalyst.cpp'; else $(CYGPATH_W) '$(srcdir)/SpectrumAnalyst.cpp'; fi`
|
||||
@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/audacity-SpectrumAnalyst.Tpo $(DEPDIR)/audacity-SpectrumAnalyst.Po
|
||||
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ $(AM_V_CXX)source='SpectrumAnalyst.cpp' object='audacity-SpectrumAnalyst.obj' libtool=no @AMDEPBACKSLASH@
|
||||
@AMDEP_TRUE@@am__fastdepCXX_FALSE@ DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
|
||||
@am__fastdepCXX_FALSE@ $(AM_V_CXX@am__nodep@)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(audacity_CPPFLAGS) $(CPPFLAGS) $(audacity_CXXFLAGS) $(CXXFLAGS) -c -o audacity-SpectrumAnalyst.obj `if test -f 'SpectrumAnalyst.cpp'; then $(CYGPATH_W) 'SpectrumAnalyst.cpp'; else $(CYGPATH_W) '$(srcdir)/SpectrumAnalyst.cpp'; fi`
|
||||
|
||||
audacity-SplashDialog.o: SplashDialog.cpp
|
||||
@am__fastdepCXX_TRUE@ $(AM_V_CXX)$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(audacity_CPPFLAGS) $(CPPFLAGS) $(audacity_CXXFLAGS) $(CXXFLAGS) -MT audacity-SplashDialog.o -MD -MP -MF $(DEPDIR)/audacity-SplashDialog.Tpo -c -o audacity-SplashDialog.o `test -f 'SplashDialog.cpp' || echo '$(srcdir)/'`SplashDialog.cpp
|
||||
@am__fastdepCXX_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/audacity-SplashDialog.Tpo $(DEPDIR)/audacity-SplashDialog.Po
|
||||
|
|
@ -9278,6 +9297,7 @@ distclean: distclean-am
|
|||
-rm -f ./$(DEPDIR)/audacity-Snap.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-SoundActivatedRecord.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-Spectrum.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-SpectrumAnalyst.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-SplashDialog.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-SseMathFuncs.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-Tags.Po
|
||||
|
|
@ -9723,6 +9743,7 @@ maintainer-clean: maintainer-clean-am
|
|||
-rm -f ./$(DEPDIR)/audacity-Snap.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-SoundActivatedRecord.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-Spectrum.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-SpectrumAnalyst.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-SplashDialog.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-SseMathFuncs.Po
|
||||
-rm -f ./$(DEPDIR)/audacity-Tags.Po
|
||||
|
|
|
|||
|
|
@ -0,0 +1,511 @@
|
|||
/**********************************************************************
|
||||
|
||||
Audacity: A Digital Audio Editor
|
||||
|
||||
SpectrumAnalyst.cpp
|
||||
|
||||
Dominic Mazzoni
|
||||
Paul Licameli split from FreqWindow.cpp
|
||||
|
||||
*******************************************************************//**
|
||||
|
||||
\class SpectrumAnalyst
|
||||
\brief Used for finding the peaks, for snapping to peaks.
|
||||
|
||||
This class is used to do the 'find peaks' snapping both in FreqPlot
|
||||
and in the spectrogram spectral selection.
|
||||
|
||||
*//*******************************************************************/
|
||||
|
||||
/*
|
||||
Salvo Ventura - November 2006
|
||||
Extended range check for additional FFT windows
|
||||
*/
|
||||
|
||||
|
||||
#include "Audacity.h"
|
||||
#include "SpectrumAnalyst.h"
|
||||
#include "FFT.h"
|
||||
|
||||
#include "SampleFormat.h"
|
||||
#include <wx/dcclient.h>
|
||||
|
||||
FreqGauge::FreqGauge(wxWindow * parent, wxWindowID winid)
|
||||
: wxStatusBar(parent, winid, wxST_SIZEGRIP)
|
||||
{
|
||||
mRange = 0;
|
||||
}
|
||||
|
||||
void FreqGauge::SetRange(int range, int bar, int gap)
|
||||
{
|
||||
mRange = range;
|
||||
mBar = bar;
|
||||
mGap = gap;
|
||||
|
||||
GetFieldRect(0, mRect);
|
||||
mRect.Inflate(-1);
|
||||
|
||||
mInterval = mRange / (mRect.width / (mBar + mGap));
|
||||
mRect.width = mBar;
|
||||
mMargin = mRect.x;
|
||||
mLast = -1;
|
||||
|
||||
Update();
|
||||
}
|
||||
|
||||
void FreqGauge::SetValue(int value)
|
||||
{
|
||||
mCur = value / mInterval;
|
||||
|
||||
if (mCur != mLast)
|
||||
{
|
||||
wxClientDC dc(this);
|
||||
dc.SetPen(*wxTRANSPARENT_PEN);
|
||||
dc.SetBrush(wxColour(100, 100, 220));
|
||||
|
||||
while (mLast < mCur)
|
||||
{
|
||||
mLast++;
|
||||
mRect.x = mMargin + mLast * (mBar + mGap);
|
||||
dc.DrawRectangle(mRect);
|
||||
}
|
||||
Update();
|
||||
}
|
||||
}
|
||||
|
||||
void FreqGauge::Reset()
|
||||
{
|
||||
mRange = 0;
|
||||
Refresh(true);
|
||||
}
|
||||
|
||||
SpectrumAnalyst::SpectrumAnalyst()
|
||||
: mAlg(Spectrum)
|
||||
, mRate(0.0)
|
||||
, mWindowSize(0)
|
||||
{
|
||||
}
|
||||
|
||||
SpectrumAnalyst::~SpectrumAnalyst()
|
||||
{
|
||||
}
|
||||
|
||||
bool SpectrumAnalyst::Calculate(Algorithm alg, int windowFunc,
|
||||
size_t windowSize, double rate,
|
||||
const float *data, size_t dataLen,
|
||||
float *pYMin, float *pYMax,
|
||||
FreqGauge *progress)
|
||||
{
|
||||
// Wipe old data
|
||||
mProcessed.resize(0);
|
||||
mRate = 0.0;
|
||||
mWindowSize = 0;
|
||||
|
||||
// Validate inputs
|
||||
int f = NumWindowFuncs();
|
||||
|
||||
if (!(windowSize >= 32 && windowSize <= 65536 &&
|
||||
alg >= SpectrumAnalyst::Spectrum &&
|
||||
alg < SpectrumAnalyst::NumAlgorithms &&
|
||||
windowFunc >= 0 && windowFunc < f)) {
|
||||
return false;
|
||||
}
|
||||
|
||||
if (dataLen < windowSize) {
|
||||
return false;
|
||||
}
|
||||
|
||||
// Now repopulate
|
||||
mRate = rate;
|
||||
mWindowSize = windowSize;
|
||||
mAlg = alg;
|
||||
|
||||
auto half = mWindowSize / 2;
|
||||
mProcessed.resize(mWindowSize);
|
||||
|
||||
Floats in{ mWindowSize };
|
||||
Floats out{ mWindowSize };
|
||||
Floats out2{ mWindowSize };
|
||||
Floats win{ mWindowSize };
|
||||
|
||||
for (size_t i = 0; i < mWindowSize; i++) {
|
||||
mProcessed[i] = 0.0f;
|
||||
win[i] = 1.0f;
|
||||
}
|
||||
|
||||
WindowFunc(windowFunc, mWindowSize, win.get());
|
||||
|
||||
// Scale window such that an amplitude of 1.0 in the time domain
|
||||
// shows an amplitude of 0dB in the frequency domain
|
||||
double wss = 0;
|
||||
for (size_t i = 0; i<mWindowSize; i++)
|
||||
wss += win[i];
|
||||
if(wss > 0)
|
||||
wss = 4.0 / (wss*wss);
|
||||
else
|
||||
wss = 1.0;
|
||||
|
||||
if (progress) {
|
||||
progress->SetRange(dataLen);
|
||||
}
|
||||
|
||||
size_t start = 0;
|
||||
int windows = 0;
|
||||
while (start + mWindowSize <= dataLen) {
|
||||
for (size_t i = 0; i < mWindowSize; i++)
|
||||
in[i] = win[i] * data[start + i];
|
||||
|
||||
switch (alg) {
|
||||
case Spectrum:
|
||||
PowerSpectrum(mWindowSize, in.get(), out.get());
|
||||
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] += out[i];
|
||||
break;
|
||||
|
||||
case Autocorrelation:
|
||||
case CubeRootAutocorrelation:
|
||||
case EnhancedAutocorrelation:
|
||||
|
||||
// Take FFT
|
||||
RealFFT(mWindowSize, in.get(), out.get(), out2.get());
|
||||
// Compute power
|
||||
for (size_t i = 0; i < mWindowSize; i++)
|
||||
in[i] = (out[i] * out[i]) + (out2[i] * out2[i]);
|
||||
|
||||
if (alg == Autocorrelation) {
|
||||
for (size_t i = 0; i < mWindowSize; i++)
|
||||
in[i] = sqrt(in[i]);
|
||||
}
|
||||
if (alg == CubeRootAutocorrelation ||
|
||||
alg == EnhancedAutocorrelation) {
|
||||
// Tolonen and Karjalainen recommend taking the cube root
|
||||
// of the power, instead of the square root
|
||||
|
||||
for (size_t i = 0; i < mWindowSize; i++)
|
||||
in[i] = pow(in[i], 1.0f / 3.0f);
|
||||
}
|
||||
// Take FFT
|
||||
RealFFT(mWindowSize, in.get(), out.get(), out2.get());
|
||||
|
||||
// Take real part of result
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] += out[i];
|
||||
break;
|
||||
|
||||
case Cepstrum:
|
||||
RealFFT(mWindowSize, in.get(), out.get(), out2.get());
|
||||
|
||||
// Compute log power
|
||||
// Set a sane lower limit assuming maximum time amplitude of 1.0
|
||||
{
|
||||
float power;
|
||||
float minpower = 1e-20*mWindowSize*mWindowSize;
|
||||
for (size_t i = 0; i < mWindowSize; i++)
|
||||
{
|
||||
power = (out[i] * out[i]) + (out2[i] * out2[i]);
|
||||
if(power < minpower)
|
||||
in[i] = log(minpower);
|
||||
else
|
||||
in[i] = log(power);
|
||||
}
|
||||
// Take IFFT
|
||||
InverseRealFFT(mWindowSize, in.get(), NULL, out.get());
|
||||
|
||||
// Take real part of result
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] += out[i];
|
||||
}
|
||||
|
||||
break;
|
||||
|
||||
default:
|
||||
wxASSERT(false);
|
||||
break;
|
||||
} //switch
|
||||
|
||||
// Update the progress bar
|
||||
if (progress) {
|
||||
progress->SetValue(start);
|
||||
}
|
||||
|
||||
start += half;
|
||||
windows++;
|
||||
}
|
||||
|
||||
if (progress) {
|
||||
// Reset for next time
|
||||
progress->Reset();
|
||||
}
|
||||
|
||||
float mYMin = 1000000, mYMax = -1000000;
|
||||
double scale;
|
||||
switch (alg) {
|
||||
case Spectrum:
|
||||
// Convert to decibels
|
||||
mYMin = 1000000.;
|
||||
mYMax = -1000000.;
|
||||
scale = wss / (double)windows;
|
||||
for (size_t i = 0; i < half; i++)
|
||||
{
|
||||
mProcessed[i] = 10 * log10(mProcessed[i] * scale);
|
||||
if(mProcessed[i] > mYMax)
|
||||
mYMax = mProcessed[i];
|
||||
else if(mProcessed[i] < mYMin)
|
||||
mYMin = mProcessed[i];
|
||||
}
|
||||
break;
|
||||
|
||||
case Autocorrelation:
|
||||
case CubeRootAutocorrelation:
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] = mProcessed[i] / windows;
|
||||
|
||||
// Find min/max
|
||||
mYMin = mProcessed[0];
|
||||
mYMax = mProcessed[0];
|
||||
for (size_t i = 1; i < half; i++)
|
||||
if (mProcessed[i] > mYMax)
|
||||
mYMax = mProcessed[i];
|
||||
else if (mProcessed[i] < mYMin)
|
||||
mYMin = mProcessed[i];
|
||||
break;
|
||||
|
||||
case EnhancedAutocorrelation:
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] = mProcessed[i] / windows;
|
||||
|
||||
// Peak Pruning as described by Tolonen and Karjalainen, 2000
|
||||
|
||||
// Clip at zero, copy to temp array
|
||||
for (size_t i = 0; i < half; i++) {
|
||||
if (mProcessed[i] < 0.0)
|
||||
mProcessed[i] = float(0.0);
|
||||
out[i] = mProcessed[i];
|
||||
}
|
||||
|
||||
// Subtract a time-doubled signal (linearly interp.) from the original
|
||||
// (clipped) signal
|
||||
for (size_t i = 0; i < half; i++)
|
||||
if ((i % 2) == 0)
|
||||
mProcessed[i] -= out[i / 2];
|
||||
else
|
||||
mProcessed[i] -= ((out[i / 2] + out[i / 2 + 1]) / 2);
|
||||
|
||||
// Clip at zero again
|
||||
for (size_t i = 0; i < half; i++)
|
||||
if (mProcessed[i] < 0.0)
|
||||
mProcessed[i] = float(0.0);
|
||||
|
||||
// Find NEW min/max
|
||||
mYMin = mProcessed[0];
|
||||
mYMax = mProcessed[0];
|
||||
for (size_t i = 1; i < half; i++)
|
||||
if (mProcessed[i] > mYMax)
|
||||
mYMax = mProcessed[i];
|
||||
else if (mProcessed[i] < mYMin)
|
||||
mYMin = mProcessed[i];
|
||||
break;
|
||||
|
||||
case Cepstrum:
|
||||
for (size_t i = 0; i < half; i++)
|
||||
mProcessed[i] = mProcessed[i] / windows;
|
||||
|
||||
// Find min/max, ignoring first and last few values
|
||||
{
|
||||
size_t ignore = 4;
|
||||
mYMin = mProcessed[ignore];
|
||||
mYMax = mProcessed[ignore];
|
||||
for (size_t i = ignore + 1; i + ignore < half; i++)
|
||||
if (mProcessed[i] > mYMax)
|
||||
mYMax = mProcessed[i];
|
||||
else if (mProcessed[i] < mYMin)
|
||||
mYMin = mProcessed[i];
|
||||
}
|
||||
break;
|
||||
|
||||
default:
|
||||
wxASSERT(false);
|
||||
break;
|
||||
}
|
||||
|
||||
if (pYMin)
|
||||
*pYMin = mYMin;
|
||||
if (pYMax)
|
||||
*pYMax = mYMax;
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
const float *SpectrumAnalyst::GetProcessed() const
|
||||
{
|
||||
return &mProcessed[0];
|
||||
}
|
||||
|
||||
int SpectrumAnalyst::GetProcessedSize() const
|
||||
{
|
||||
return mProcessed.size() / 2;
|
||||
}
|
||||
|
||||
float SpectrumAnalyst::GetProcessedValue(float freq0, float freq1) const
|
||||
{
|
||||
float bin0, bin1, binwidth;
|
||||
|
||||
if (mAlg == Spectrum) {
|
||||
bin0 = freq0 * mWindowSize / mRate;
|
||||
bin1 = freq1 * mWindowSize / mRate;
|
||||
} else {
|
||||
bin0 = freq0 * mRate;
|
||||
bin1 = freq1 * mRate;
|
||||
}
|
||||
binwidth = bin1 - bin0;
|
||||
|
||||
float value = float(0.0);
|
||||
|
||||
if (binwidth < 1.0) {
|
||||
float binmid = (bin0 + bin1) / 2.0;
|
||||
int ibin = (int)(binmid) - 1;
|
||||
if (ibin < 1)
|
||||
ibin = 1;
|
||||
if (ibin >= GetProcessedSize() - 3)
|
||||
ibin = std::max(0, GetProcessedSize() - 4);
|
||||
|
||||
value = CubicInterpolate(mProcessed[ibin],
|
||||
mProcessed[ibin + 1],
|
||||
mProcessed[ibin + 2],
|
||||
mProcessed[ibin + 3], binmid - ibin);
|
||||
|
||||
} else {
|
||||
if (bin0 < 0)
|
||||
bin0 = 0;
|
||||
if (bin1 >= GetProcessedSize())
|
||||
bin1 = GetProcessedSize() - 1;
|
||||
|
||||
if ((int)(bin1) > (int)(bin0))
|
||||
value += mProcessed[(int)(bin0)] * ((int)(bin0) + 1 - bin0);
|
||||
bin0 = 1 + (int)(bin0);
|
||||
while (bin0 < (int)(bin1)) {
|
||||
value += mProcessed[(int)(bin0)];
|
||||
bin0 += 1.0;
|
||||
}
|
||||
value += mProcessed[(int)(bin1)] * (bin1 - (int)(bin1));
|
||||
|
||||
value /= binwidth;
|
||||
}
|
||||
|
||||
return value;
|
||||
}
|
||||
|
||||
float SpectrumAnalyst::FindPeak(float xPos, float *pY) const
|
||||
{
|
||||
float bestpeak = 0.0f;
|
||||
float bestValue = 0.0;
|
||||
if (GetProcessedSize() > 1) {
|
||||
bool up = (mProcessed[1] > mProcessed[0]);
|
||||
float bestdist = 1000000;
|
||||
for (int bin = 3; bin < GetProcessedSize() - 1; bin++) {
|
||||
bool nowUp = mProcessed[bin] > mProcessed[bin - 1];
|
||||
if (!nowUp && up) {
|
||||
// Local maximum. Find actual value by cubic interpolation
|
||||
int leftbin = bin - 2;
|
||||
/*
|
||||
if (leftbin < 1)
|
||||
leftbin = 1;
|
||||
*/
|
||||
float valueAtMax = 0.0;
|
||||
float max = leftbin + CubicMaximize(mProcessed[leftbin],
|
||||
mProcessed[leftbin + 1],
|
||||
mProcessed[leftbin + 2],
|
||||
mProcessed[leftbin + 3],
|
||||
&valueAtMax);
|
||||
|
||||
float thispeak;
|
||||
if (mAlg == Spectrum)
|
||||
thispeak = max * mRate / mWindowSize;
|
||||
else
|
||||
thispeak = max / mRate;
|
||||
|
||||
if (fabs(thispeak - xPos) < bestdist) {
|
||||
bestpeak = thispeak;
|
||||
bestdist = fabs(thispeak - xPos);
|
||||
bestValue = valueAtMax;
|
||||
// Should this test come after the enclosing if?
|
||||
if (thispeak > xPos)
|
||||
break;
|
||||
}
|
||||
}
|
||||
up = nowUp;
|
||||
}
|
||||
}
|
||||
|
||||
if (pY)
|
||||
*pY = bestValue;
|
||||
return bestpeak;
|
||||
}
|
||||
|
||||
// If f(0)=y0, f(1)=y1, f(2)=y2, and f(3)=y3, this function finds
|
||||
// the degree-three polynomial which best fits these points and
|
||||
// returns the value of this polynomial at a value x. Usually
|
||||
// 0 < x < 3
|
||||
|
||||
float SpectrumAnalyst::CubicInterpolate(float y0, float y1, float y2, float y3, float x) const
|
||||
{
|
||||
float a, b, c, d;
|
||||
|
||||
a = y0 / -6.0 + y1 / 2.0 - y2 / 2.0 + y3 / 6.0;
|
||||
b = y0 - 5.0 * y1 / 2.0 + 2.0 * y2 - y3 / 2.0;
|
||||
c = -11.0 * y0 / 6.0 + 3.0 * y1 - 3.0 * y2 / 2.0 + y3 / 3.0;
|
||||
d = y0;
|
||||
|
||||
float xx = x * x;
|
||||
float xxx = xx * x;
|
||||
|
||||
return (a * xxx + b * xx + c * x + d);
|
||||
}
|
||||
|
||||
float SpectrumAnalyst::CubicMaximize(float y0, float y1, float y2, float y3, float * max) const
|
||||
{
|
||||
// Find coefficients of cubic
|
||||
|
||||
float a, b, c, d;
|
||||
|
||||
a = y0 / -6.0 + y1 / 2.0 - y2 / 2.0 + y3 / 6.0;
|
||||
b = y0 - 5.0 * y1 / 2.0 + 2.0 * y2 - y3 / 2.0;
|
||||
c = -11.0 * y0 / 6.0 + 3.0 * y1 - 3.0 * y2 / 2.0 + y3 / 3.0;
|
||||
d = y0;
|
||||
|
||||
// Take derivative
|
||||
|
||||
float da, db, dc;
|
||||
|
||||
da = 3 * a;
|
||||
db = 2 * b;
|
||||
dc = c;
|
||||
|
||||
// Find zeroes of derivative using quadratic equation
|
||||
|
||||
float discriminant = db * db - 4 * da * dc;
|
||||
if (discriminant < 0.0)
|
||||
return float(-1.0); // error
|
||||
|
||||
float x1 = (-db + sqrt(discriminant)) / (2 * da);
|
||||
float x2 = (-db - sqrt(discriminant)) / (2 * da);
|
||||
|
||||
// The one which corresponds to a local _maximum_ in the
|
||||
// cubic is the one we want - the one with a negative
|
||||
// second derivative
|
||||
|
||||
float dda = 2 * da;
|
||||
float ddb = db;
|
||||
|
||||
if (dda * x1 + ddb < 0)
|
||||
{
|
||||
*max = a*x1*x1*x1+b*x1*x1+c*x1+d;
|
||||
return x1;
|
||||
}
|
||||
else
|
||||
{
|
||||
*max = a*x2*x2*x2+b*x2*x2+c*x2+d;
|
||||
return x2;
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,82 @@
|
|||
/**********************************************************************
|
||||
|
||||
Audacity: A Digital Audio Editor
|
||||
|
||||
SpectrumAnalyst.h
|
||||
|
||||
Dominic Mazzoni
|
||||
Paul Licameli split from FreqWindow.h
|
||||
|
||||
**********************************************************************/
|
||||
|
||||
#ifndef __AUDACITY_SPECTRUM_ANALYST__
|
||||
#define __AUDACITY_SPECTRUM_ANALYST__
|
||||
|
||||
#include <vector>
|
||||
#include <wx/statusbr.h>
|
||||
|
||||
class FreqGauge;
|
||||
|
||||
class SpectrumAnalyst
|
||||
{
|
||||
public:
|
||||
|
||||
enum Algorithm {
|
||||
Spectrum,
|
||||
Autocorrelation,
|
||||
CubeRootAutocorrelation,
|
||||
EnhancedAutocorrelation,
|
||||
Cepstrum,
|
||||
|
||||
NumAlgorithms
|
||||
};
|
||||
|
||||
SpectrumAnalyst();
|
||||
~SpectrumAnalyst();
|
||||
|
||||
// Return true iff successful
|
||||
bool Calculate(Algorithm alg,
|
||||
int windowFunc, // see FFT.h for values
|
||||
size_t windowSize, double rate,
|
||||
const float *data, size_t dataLen,
|
||||
float *pYMin = NULL, float *pYMax = NULL, // outputs
|
||||
FreqGauge *progress = NULL);
|
||||
|
||||
const float *GetProcessed() const;
|
||||
int GetProcessedSize() const;
|
||||
|
||||
float GetProcessedValue(float freq0, float freq1) const;
|
||||
float FindPeak(float xPos, float *pY) const;
|
||||
|
||||
private:
|
||||
float CubicInterpolate(float y0, float y1, float y2, float y3, float x) const;
|
||||
float CubicMaximize(float y0, float y1, float y2, float y3, float * max) const;
|
||||
|
||||
private:
|
||||
Algorithm mAlg;
|
||||
double mRate;
|
||||
size_t mWindowSize;
|
||||
std::vector<float> mProcessed;
|
||||
};
|
||||
|
||||
class FreqGauge final : public wxStatusBar
|
||||
{
|
||||
public:
|
||||
FreqGauge(wxWindow * parent, wxWindowID winid);
|
||||
|
||||
void SetRange(int range, int bar = 12, int gap = 3);
|
||||
void SetValue(int value);
|
||||
void Reset();
|
||||
|
||||
private:
|
||||
wxRect mRect;
|
||||
int mRange;
|
||||
int mCur;
|
||||
int mLast;
|
||||
int mInterval;
|
||||
int mBar;
|
||||
int mGap;
|
||||
int mMargin;
|
||||
};
|
||||
|
||||
#endif
|
||||
|
|
@ -4,7 +4,7 @@
|
|||
#include "../AdornedRulerPanel.h"
|
||||
#include "../AudioIO.h"
|
||||
#include "../CommonCommandFlags.h"
|
||||
#include "../FreqWindow.h"
|
||||
#include "../SpectrumAnalyst.h"
|
||||
#include "../Prefs.h"
|
||||
#include "../Project.h"
|
||||
#include "../ProjectAudioIO.h"
|
||||
|
|
|
|||
|
|
@ -17,7 +17,7 @@ Paul Licameli split from TrackPanel.cpp
|
|||
#include "TrackView.h"
|
||||
|
||||
#include "../../AColor.h"
|
||||
#include "../../FreqWindow.h"
|
||||
#include "../../SpectrumAnalyst.h"
|
||||
#include "../../NumberScale.h"
|
||||
#include "../../Project.h"
|
||||
#include "../../ProjectAudioIO.h"
|
||||
|
|
|
|||
Loading…
Reference in New Issue