Skip to content

Commit

Permalink
Merge pull request #50 from arthurits/main
Browse files Browse the repository at this point in the history
FFTfreq: fix off-by-one error
  • Loading branch information
swharden committed Jun 16, 2022
2 parents 454708e + 0809ad7 commit 77c56f3
Show file tree
Hide file tree
Showing 8 changed files with 164 additions and 34 deletions.
Binary file modified dev/quickstart/fft-windowed.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified dev/quickstart/fft.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified dev/quickstart/periodogram.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
78 changes: 78 additions & 0 deletions src/FftSharp.Tests/FftFreqTests.cs
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
using NUnit.Framework;
using System;
using System.Linq;

namespace FftSharp.Tests
{
internal class FftFreqTests
{
[Test]
public void Test_FftFreq_KnownValues()
{
// https://github.com/swharden/FftSharp/issues/49

// generate signal with 2 Hz sine wave
int sampleCount = 16;
double sampleRateHz = 16;
double sinFrequencyHz = 2;
double[] samples = Enumerable.Range(0, sampleCount)
.Select(i => Math.Sin(2 * Math.PI * sinFrequencyHz * i / sampleRateHz))
.ToArray();

// get FFT
double[] fft = FftSharp.Transform.FFTmagnitude(samples);
double[] fftKnown = { 0, 0, 1, 0, 0, 0, 0, 0, 0 };
Assert.That(fft, Is.EqualTo(fftKnown).Within(1e-10));

// calculate FFT frequencies both ways
double[] fftFreqKnown = { 0, 1, 2, 3, 4, 5, 6, 7, 8 };
double[] fftFreq = FftSharp.Transform.FFTfreq(sampleRateHz, fft.Length);
double[] fftFreq2 = FftSharp.Transform.FFTfreq(sampleRateHz, fft);
Assert.That(fftFreq, Is.EqualTo(fftFreqKnown));
Assert.That(fftFreq2, Is.EqualTo(fftFreqKnown));

ScottPlot.Plot plt1 = new(400, 200);
plt1.AddSignal(samples, sampleRateHz);
TestTools.SaveFig(plt1, "signal");

ScottPlot.Plot plt2 = new(400, 200);
plt2.AddScatter(fftFreq, fft);
plt2.AddVerticalLine(2, System.Drawing.Color.Red, style: ScottPlot.LineStyle.Dash);
TestTools.SaveFig(plt2, "fft");
}

[Test]
public void Test_Freq_Lookup()
{
/* Compare against values generated using Python:
*
* >>> numpy.fft.fftfreq(16, 1/16000)
* array([ 0., 1000., 2000., 3000., 4000., 5000., 6000., 7000.,
* -8000., -7000., -6000., -5000., -4000., -3000., -2000., -1000.])
*
*/

int sampleRate = 16000;
int pointCount = 16;
double[] signal = new double[pointCount];
double[] fft = FftSharp.Transform.FFTmagnitude(signal);

double[] freqsFullKnown = {
0, 1000, 2000, 3000, 4000, 5000, 6000, 7000,
-8000, -7000, -6000, -5000, -4000, -3000, -2000, -1000
};
double[] freqsFull = FftSharp.Transform.FFTfreq(sampleRate, signal, oneSided: false);
double[] freqsFull2 = FftSharp.Transform.FFTfreq(sampleRate, signal.Length, oneSided: false);
Assert.That(freqsFull, Is.EqualTo(freqsFullKnown));
Assert.That(freqsFull2, Is.EqualTo(freqsFullKnown));

double[] freqsHalfKnown = {
0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000
};
double[] freqsHalf = FftSharp.Transform.FFTfreq(sampleRate, fft, oneSided: true);
double[] freqsHalf2 = FftSharp.Transform.FFTfreq(sampleRate, fft.Length, oneSided: true);
Assert.That(freqsHalf, Is.EqualTo(freqsHalfKnown));
Assert.That(freqsHalf2, Is.EqualTo(freqsHalfKnown));
}
}
}
28 changes: 0 additions & 28 deletions src/FftSharp.Tests/FreqLookup.cs

This file was deleted.

63 changes: 63 additions & 0 deletions src/FftSharp.Tests/Transform.cs
Original file line number Diff line number Diff line change
Expand Up @@ -155,5 +155,68 @@ public void Test_FftInput_ThrowsIfLengthIsNotPowerOfTwo(int length, bool shouldT
Assert.DoesNotThrow(realSpanRFFT);
}
}

[TestCase(0, true)]
[TestCase(1, false)]
[TestCase(2, false)]
[TestCase(4, false)]
[TestCase(8, false)]
[TestCase(16, false)]
[TestCase(32, false)]
[TestCase(64, false)]
public void Test_Fft_Complex_DifferentLengths(int pointCount, bool shouldFail)
{
Complex[] signal = new Complex[pointCount];
if (shouldFail)
{
Assert.Throws<ArgumentException>(() => FftSharp.Transform.FFT(signal));
}
else
{
Assert.DoesNotThrow(() => FftSharp.Transform.FFT(signal));
}
}

[TestCase(0, true)]
[TestCase(1, false)]
[TestCase(2, false)]
[TestCase(4, false)]
[TestCase(8, false)]
[TestCase(16, false)]
[TestCase(32, false)]
[TestCase(64, false)]
public void Test_Fft_Double_DifferentLengths(int pointCount, bool shouldFail)
{
double[] signal = new double[pointCount];
if (shouldFail)
{
Assert.Throws<ArgumentException>(() => FftSharp.Transform.FFT(signal));
}
else
{
Assert.DoesNotThrow(() => FftSharp.Transform.FFT(signal));
}
}

[TestCase(0, true)]
[TestCase(1, true)]
[TestCase(2, true)]
[TestCase(4, true)]
[TestCase(8, true)]
[TestCase(16, false)]
[TestCase(32, false)]
[TestCase(64, false)]
public void Test_Fft_Magnitude_DifferentLengths(int pointCount, bool shouldFail)
{
double[] signal = new double[pointCount];
if (shouldFail)
{
Assert.Throws<ArgumentException>(() => FftSharp.Transform.FFTmagnitude(signal));
}
else
{
Assert.DoesNotThrow(() => FftSharp.Transform.FFTmagnitude(signal));
}
}
}
}
10 changes: 5 additions & 5 deletions src/FftSharp/FftSharp.csproj
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,9 @@
<PackageReleaseNotes>https://github.com/swharden/FftSharp/releases</PackageReleaseNotes>
<IncludeSymbols>true</IncludeSymbols>
<DocumentationFile>FftSharp.xml</DocumentationFile>
<Version>1.1.5</Version>
<AssemblyVersion>1.1.5.0</AssemblyVersion>
<FileVersion>1.1.5.0</FileVersion>
<Version>1.1.6</Version>
<AssemblyVersion>1.1.6.0</AssemblyVersion>
<FileVersion>1.1.6.0</FileVersion>
<NoWarn>1591</NoWarn>
<DebugType>portable</DebugType>
<IncludeSymbols>true</IncludeSymbols>
Expand All @@ -34,15 +34,15 @@
</ItemGroup>

<ItemGroup>
<PackageReference Include="Microsoft.SourceLink.GitHub" Version="1.0.0">
<PackageReference Include="Microsoft.SourceLink.GitHub" Version="1.1.1">
<IncludeAssets>runtime; build; native; contentfiles; analyzers; buildtransitive</IncludeAssets>
<PrivateAssets>all</PrivateAssets>
</PackageReference>
</ItemGroup>

<ItemGroup Condition="'$(TargetFramework)' == 'netstandard2.0'">
<PackageReference Include="System.ValueTuple" Version="4.5.0" />
<PackageReference Include="System.Memory" Version="4.5.4" />
<PackageReference Include="System.Memory" Version="4.5.5" />
</ItemGroup>

</Project>
19 changes: 18 additions & 1 deletion src/FftSharp/Transform.cs
Original file line number Diff line number Diff line change
Expand Up @@ -122,13 +122,16 @@ private static int BitReverse(int value, int maxValue)
/// <summary>
/// Calculate sample frequency for each point in a FFT
/// </summary>
/// <param name="sampleRate">Sample rate (Hz) of the original signal</param>
/// <param name="pointCount">Number of points to generate (typically the length of the FFT)</param>
/// <param name="oneSided">Whether or not frequencies are for a one-sided FFT (containing only real numbers)</param>
public static double[] FFTfreq(double sampleRate, int pointCount, bool oneSided = true)
{
double[] freqs = new double[pointCount];

if (oneSided)
{
double fftPeriodHz = sampleRate / pointCount / 2;
double fftPeriodHz = sampleRate / (pointCount - 1) / 2;

// freqs start at 0 and approach maxFreq
for (int i = 0; i < pointCount; i++)
Expand All @@ -151,6 +154,17 @@ public static double[] FFTfreq(double sampleRate, int pointCount, bool oneSided
}
}

/// <summary>
/// Calculate sample frequency for each point in a FFT
/// </summary>
/// <param name="sampleRate">Sample rate (Hz) of the original signal</param>
/// <param name="fft">FFT array for which frequencies should be generated</param>
/// <param name="oneSided">Whether or not frequencies are for a one-sided FFT (containing only real numbers)</param>
public static double[] FFTfreq(double sampleRate, double[] fft, bool oneSided = true)
{
return FFTfreq(sampleRate, fft.Length, oneSided);
}

/// <summary>
/// Return the distance between each FFT point in frequency units (Hz)
/// </summary>
Expand Down Expand Up @@ -294,6 +308,9 @@ public static double[] FFTmagnitude(double[] input)
/// <param name="input">real input</param>
public static void FFTmagnitude(Span<double> destination, Span<double> input)
{
if (input.Length < 16)
throw new ArgumentException("This overload requires an input with at least 16 points");

if (!IsPowerOfTwo(input.Length))
throw new ArgumentException("Input length must be an even power of 2");

Expand Down

0 comments on commit 77c56f3

Please sign in to comment.