Attempting to replicate Simulink's Fourier Block for variable sample times in MATLAB.

I am working on replicating the Simulink Fourier Block's functionality within MATLAB for real-time signal analysis, focusing on a fundamental frequency of 50Hz. The challenge arises when handling varying sample times, e.g., 1e-05 vs. 1.5e-05 seconds. While the original Simulink block functions correctly under any sample time, my MATLAB implementation struggles with non-integer numbers in the buffer length, crucial for signal discretization. Attempts to adjust, such as dynamic buffer length and windowing techniques, have been made but fail to provide a universal solution. I'm looking for advice or methods to overcome this challenge and ensure accurate Fourier analysis for any sample time. Here's the code I'm working with:

function [magnitude, phase] = RealTimeFourier(inputSignal, fundamentalFreq, harmonicN, sampleTime, initialMagnitude, initialPhase)

%#codegen

persistent signalBuffer cosBuffer sinBuffer bufferIndex hasCompletedCycle

if isempty(signalBuffer)

bufferSize = ceil(1 / (fundamentalFreq * sampleTime));

signalBuffer = zeros(1, bufferSize);

cosBuffer = cos(2 * pi * fundamentalFreq * harmonicN * (0:bufferSize-1) * sampleTime);

sinBuffer = sin(2 * pi * fundamentalFreq * harmonicN * (0:bufferSize-1) * sampleTime);

bufferIndex = 1;

hasCompletedCycle = false;

end

signalBuffer(bufferIndex) = inputSignal;

if bufferIndex == length(signalBuffer) && ~hasCompletedCycle

hasCompletedCycle = true;

end

bufferIndex = mod(bufferIndex, length(signalBuffer)) + 1;

if hasCompletedCycle

a_n = 2 / length(signalBuffer) * sum(signalBuffer .* cosBuffer);

b_n = 2 / length(signalBuffer) * sum(signalBuffer .* sinBuffer);

magnitude = sqrt(a_n^2 + b_n^2);

phase = atand(a_n / b_n);

else

magnitude = initialMagnitude;

phase = initialPhase;

end

end

Answer 1

[Namnendra]

Hi Benjamin,

To handle varying sample times accurately in your real-time Fourier analysis, you need to dynamically adjust the buffer size and ensure that your calculations remain accurate regardless of the sample time. Here are a few suggestions to improve your implementation:

1. Dynamic Buffer Size Adjustment: Calculate the buffer size based on the sample time and fundamental frequency to ensure it covers exactly one period of the fundamental frequency.

2. Windowing Technique: Apply a windowing technique to reduce spectral leakage, which can be significant when the buffer size is not an integer multiple of the signal period.

3. Accurate Phase Calculation: Ensure the phase calculation handles the quadrant correctly by using `atan2` instead of `atand`.

Here's a revised version of your function incorporating these suggestions:

function [magnitude, phase] = RealTimeFourier(inputSignal, fundamentalFreq, harmonicN, sampleTime, initialMagnitude, initialPhase)
%#codegen

persistent signalBuffer cosBuffer sinBuffer bufferSize bufferIndex hasCompletedCycle
if isempty(signalBuffer)
    % Calculate buffer size to cover exactly one period of the fundamental frequency
    bufferSize = round(1 / (fundamentalFreq * sampleTime));
    signalBuffer = zeros(1, bufferSize);
    cosBuffer = cos(2 * pi * harmonicN * (0:bufferSize-1) * sampleTime);
    sinBuffer = sin(2 * pi * harmonicN * (0:bufferSize-1) * sampleTime);
    bufferIndex = 1;
    hasCompletedCycle = false;
end

% Update the signal buffer
signalBuffer(bufferIndex) = inputSignal;

% Check if the buffer has completed a cycle
if bufferIndex == bufferSize
    hasCompletedCycle = true;
end

% Update the buffer index
bufferIndex = mod(bufferIndex, bufferSize) + 1;

if hasCompletedCycle
    % Apply a windowing function (Hann window) to reduce spectral leakage
    window = hann(bufferSize)';
    windowedSignal = signalBuffer .* window;
    
    % Perform the Fourier analysis
    a_n = 2 / bufferSize * sum(windowedSignal .* cosBuffer);
    b_n = 2 / bufferSize * sum(windowedSignal .* sinBuffer);
    magnitude = sqrt(a_n^2 + b_n^2);
    phase = atan2(b_n, a_n);  % Use atan2 to handle the quadrant correctly
else
    magnitude = initialMagnitude;
    phase = initialPhase;
end

end

Explanation

1. Dynamic Buffer Size Calculation: The buffer size is calculated using `round` to ensure it covers exactly one period of the fundamental frequency, regardless of the sample time.

2. Windowing Function: A Hann window is applied to the signal buffer to reduce spectral leakage. This is particularly useful when the buffer size is not an integer multiple of the signal period.

3. Accurate Phase Calculation: The phase is calculated using `atan2(b_n, a_n)` to handle the correct quadrant of the phase angle.

Usage

You can use this function in a loop to process your real-time signal data:

% Example usage
fs = 1e5; % Sample rate
T = 1/fs; % Sample time
f0 = 50; % Fundamental frequency
harmonicN = 1; % Harmonic number
initialMagnitude = 0;
initialPhase = 0;

% Simulate a signal
t = 0:T:1; % Time vector
inputSignal = sin(2 * pi * f0 * t);

% Initialize output arrays
magnitude = zeros(size(t));
phase = zeros(size(t));

% Process the signal in real-time
for i = 1:length(t)
    [magnitude(i), phase(i)] = RealTimeFourier(inputSignal(i), f0, harmonicN, T, initialMagnitude, initialPhase);
end

% Plot the results
figure;
subplot(2, 1, 1);
plot(t, magnitude);
title('Magnitude');
subplot(2, 1, 2);
plot(t, phase);
title('Phase');

This revised function should handle varying sample times more robustly and provide accurate Fourier analysis for your real-time signal processing needs.

Thank you.