swharden/FftSharp

Welch's method: mean FFT of overlapping windowed segments

Closed this issue · 10 comments

busranlu commented

I want to ask about Welch's method for better results(based on windowing with hamming etc. and sliding window with overlap)
With FftSharp, Can I apply this method? Can I adjust the window length, overlap ratio?

Thank you

swharden commented

Hi @busranlu, thanks for this question! FftSharp has window classes (including Welsh) for creating a Welch window and applying it to a segment of data. Here's some code you may find helpful:

Generating a Welch Window

FftSharp.Windows.Welch window = new();
double[] windowYs = window.Create(100);

ScottPlot.Plot plot = new();
plot.Add.Signal(windowYs);
plot.SavePng("welch.png", 600, 400).LaunchFile();

welch

Applying a Welch Window to Data

FftSharp.Windows.Welch window = new();
double[] data = ScottPlot.Generate.RandomSample(100);
double[] windowedData = window.Apply(data);

ScottPlot.Plot plot = new();
plot.Add.Signal(data);
plot.SavePng("data.png", 600, 400).LaunchFile();

plot.Clear();
plot.Add.Signal(windowedData);
plot.SavePng("windowed.png", 600, 400).LaunchFile();
data windowed
data windowed

However, I don't think FftSharp has built-in tooling for managing "sliding window" analysis of continuous signals. Some of this functionality is available in my Spectrogram library, and you may find some of the code examples there helpful even if your target application isn't creating spectrograms. https://github.com/swharden/Spectrogram

Let me know if I misunderstood your question! Good luck 👍 🚀

busranlu commented

Hi Scott,
Thanks for your answer. Maybe, I should mention the MATLAB function "pwelch". In this function, window length, overlap(50%, 25%...), nfft(64,128...), and window type(hanning, hamming..) are defined for Welch's overlapped segment averaging estimator. I wanted to ask if it is possible or not.
thank you

swharden commented

Hi @busranlu, thanks for following-up!

At first glance it looks like MATLAB's pwelch function https://www.mathworks.com/help/signal/ref/pwelch.html takes a signal, applies the Welch window, then gets the FFT of it. It has overloads that allow it to automatically "slide across a signal", doing that multiple times.

FftSharp doesn't have that sliding functionality built in, but you could probably achieve a similar affect with a for loop. I'll try to generate a code sample that mimics the forum post https://www.mathworks.com/matlabcentral/answers/1876767-how-to-use-pwelch-for-psd which may be the type of analysis you're looking for. Let me know if I'm misunderstanding your goal though!

image

busranlu commented

Dear Scott,
I agree, the loop can solve the problem. Thank you for your time and attention. This image may explain the purpose better.
image
https://www.techscience.com/cmc/v67n3/41590/html#:~:text=The%20Welch%20method%20is%20carried,of%20the%20individual%20power%20measurements.

swharden commented

Hi @busranlu, thanks for the clarification! Here's some code which may be helpful to you or others finding this discussion later. I'm reducing the overlap to only 20% to make it more obvious what's happening:

Original Signal

plot1

Moving Windows

plot2

Moving Window Segments

plot3

Moving Window Filtered Signal

plot4

Code

// show sample data
int dataSize = 500;
double[] data = ScottPlot.Generate.RandomNormal(dataSize);

ScottPlot.Plot plot1 = new();
plot1.Add.Signal(data);
plot1.Axes.SetLimitsY(-4, 4);
plot1.SavePng("plot1.png", 800, 300).LaunchFile();

// show sliding windows
int windowSize = 100;
double overlapFraction = 0.2;
int windowStep = (int)(windowSize * (1 - overlapFraction));
FftSharp.Windows.Welch window = new();
double[] windowYs = window.Create(windowSize);
ScottPlot.Plot plot2 = new();
for (int i = 0; i < dataSize - windowSize; i += windowStep)
{
    var sig = plot2.Add.Signal(windowYs);
    sig.Data.XOffset = i;
    sig.LineWidth = 2;
}
plot2.SavePng("plot2.png", 800, 300).LaunchFile();

// show sliding windows applied to data segments
double[] summedData = new double[dataSize];
ScottPlot.Plot plot3 = new();
for (int i = 0; i < dataSize - windowSize; i += windowStep)
{
    // apply the window to a segment
    double[] segment = data[i..(i + windowSize)];
    window.ApplyInPlace(segment);

    // show the windowed segment
    var sig = plot3.Add.Signal(segment);
    sig.Data.XOffset = i;
    sig.LineWidth = 2;

    // add that segment to our summed signal
    for (int j = 0; j < segment.Length; j++)
        summedData[i + j] += segment[j];
}
plot3.Axes.SetLimitsY(-4, 4);
plot3.SavePng("plot3.png", 800, 300).LaunchFile();


// show summed windowed data
ScottPlot.Plot plot4 = new();
plot4.Add.Signal(summedData);
plot4.Axes.SetLimitsY(-4, 4);
plot4.SavePng("plot4.png", 800, 300).LaunchFile();
busranlu commented

I think this code is suitable for filtering but not for obtaining the power spectrum of a signal because we are still in the time domain, not the frequency domain.

swharden commented

I think this code is suitable for filtering but not for obtaining the power spectrum of a signal because we are still in the time domain, not the frequency domain.

Good point! I'll modify the code sample above as shown below to demonstrate this strategy.

// show sliding window FFTs
ScottPlot.Plot plot3b = new();
for (int i = 0; i < dataSize - windowSize; i += windowStep)
{
    // apply the window to a segment
    double[] segment = data[i..(i + windowSize)];
    window.ApplyInPlace(segment);
    double[] segment2 = FftSharp.Pad.ZeroPad(segment);
    System.Numerics.Complex[] fft = FftSharp.FFT.Forward(segment2);
    double[] power = FftSharp.FFT.Power(fft);

    // show the windowed segment
    var sig = plot3b.Add.Signal(power);
    sig.LineWidth = 2;

    // add that segment to our summed signal
    for (int j = 0; j < segment.Length; j++)
        summedData[i + j] += segment[j];
}
plot3b.SavePng("plot3b.png", 800, 300).LaunchFile();

the PSD isn't too interesting because the source data is random noise

plot3b

busranlu commented

Dear Scott,
Firstly I want to thank you for your time and attention. This code gives us the power spectrum of every segment separately. Pwelch works differently and converts one signal as a power spectrum output despite the segmentation process. I think I understand the process and general idea.
Thank you for your help

swharden commented

Hi @busranlu, I'm glad you found the code examples helpful!

This code gives us the power spectrum of every segment separately. Pwelch works differently and converts one signal as a power spectrum output despite the segmentation process.

From the paper: The Welch method is carried out by dividing the time signal into successive overlapping segments, forming the periodogram for each block, and averaging

The code example above would be complete if we averaged the final signals together, rather than showing them as overlapping signals. If we loop over the individual FFTs and creating a new average curve, this is the result:

plot5

Full source code 
// show sample data
int dataSize = 500;
double[] data = ScottPlot.Generate.RandomNormal(dataSize);

ScottPlot.Plot plot1 = new();
plot1.Add.Signal(data);
plot1.Axes.SetLimitsY(-4, 4);
plot1.SavePng("plot1.png", 800, 300).LaunchFile();
plot1.Title("Original Data");

// show sliding windows
int windowSize = 100;
double overlapFraction = 0.2;
int windowStep = (int)(windowSize * (1 - overlapFraction));
FftSharp.Windows.Welch window = new();
double[] windowYs = window.Create(windowSize);
ScottPlot.Plot plot2 = new();
for (int i = 0; i < dataSize - windowSize; i += windowStep)
{
    var sig = plot2.Add.Signal(windowYs);
    sig.Data.XOffset = i;
    sig.LineWidth = 2;
}
plot2.Title("Segment windows");
plot2.SavePng("plot2.png", 800, 300).LaunchFile();

// show sliding windows applied to data segments
double[] summedData = new double[dataSize];
ScottPlot.Plot plot3 = new();
for (int i = 0; i < dataSize - windowSize; i += windowStep)
{
    // apply the window to a segment
    double[] segment = data[i..(i + windowSize)];
    window.ApplyInPlace(segment);

    // show the windowed segment
    var sig = plot3.Add.Signal(segment);
    sig.Data.XOffset = i;
    sig.LineWidth = 2;

    // add that segment to our summed signal
    for (int j = 0; j < segment.Length; j++)
        summedData[i + j] += segment[j];
}
plot3.Axes.SetLimitsY(-4, 4);
plot3.Title("Individual windowed segments");
plot3.SavePng("plot3.png", 800, 300).LaunchFile();

// show sliding window FFTs
ScottPlot.Plot plot3b = new();
List<double[]> ffts = [];
for (int i = 0; i < dataSize - windowSize; i += windowStep)
{
    // apply the window to a segment
    double[] segment = data[i..(i + windowSize)];
    window.ApplyInPlace(segment);
    double[] segment2 = FftSharp.Pad.ZeroPad(segment);
    System.Numerics.Complex[] fft = FftSharp.FFT.Forward(segment2);
    double[] power = FftSharp.FFT.Power(fft);
    ffts.Add(power);

    // show the windowed segment
    var sig = plot3b.Add.Signal(power);
    sig.LineWidth = 2;

    // add that segment to our summed signal
    for (int j = 0; j < segment.Length; j++)
        summedData[i + j] += segment[j];
}
plot3b.Title("FFT of individual windowed segments");
plot3b.SavePng("plot3b.png", 800, 300).LaunchFile();

// show summed windowed data
ScottPlot.Plot plot4 = new();
plot4.Add.Signal(summedData);
plot4.Axes.SetLimitsY(-4, 4);
plot4.Title("Sum of all windowed segments");
plot4.SavePng("plot4.png", 800, 300).LaunchFile();

// show summed windowed data
double[] meanFft = new double[ffts.First().Length];
for (int i = 0; i < ffts.Count; i++)
{
    for (int j = 0; j < ffts[i].Length; j++)
    {
        meanFft[j] += ffts[i][j] / ffts.Count;
    }
}
ScottPlot.Plot plot5 = new();
plot5.Add.Signal(meanFft);
plot5.Title("Mean FFT for all segments");
plot5.SavePng("plot5.png", 800, 300).LaunchFile();
swharden commented

here's the code from the previous comment applied to a signal containing a noisy sine wave

plot1

plot3

plot3b

plot5