Describing my graph in terms of a fourier series.

So I have a graph as shown in the attached image. It represents the transient velocity of blood flow during one cardiac cycle. it repeats like this. I am needing to write some code (to input this velocity function into my cfd ANSYS model) and i want to obtain a Fourier series function that will describe the graph in question. How can I do that automatically on matlab? I used the grabit extension to draw this graph from an image I have - so I have all the data points of the curve.

Thanks!

[Les Beckham]

That data has a lot of discontinuities in the first derivative which will make this problem quite difficult. If, as you say, the data simply repeats (for many cycles?), you might try using a spline fit over several cycles. However, you have not shown even one complete cycle as your start and end values are very different.

You could also try smoothing the data (moving average, savitzky-golay?). It looks like relatively low sample rate relative to the higher frequencies in the data. How much of the higher frequency content do you need to retain?

Why do you need a fourier series? That may be the most pertinent question here. Sometimes, I notice that if I assume an approach to a problem initially, I get frustrated with making that approach work when, if I hadn't made that assumption in the first place and explored other options initially, I might have found the right approach sooner.

Look again at your data over a longer time period that actually shows the long term periodic behavior (I hope you have a more complete data set) and keep your mind open to other approaches than a fourier series.

So, after that ramble, please respond with answers to these questions and we might be able to help:

1. Why do you think that you need a fourier series to fit the data?
2. Can we really assume that the data is periodic over a longer time period?
3. Do you have a longer dataset that shows the longer term periodic behavior?

If the answer to #3 above is yes, please attach the data to your original question. That would help a lot.

[Walter Roberson]

The start and end values look decently similar to me, both roughly 0.1?

[Les Beckham]

I stand corrected. Thanks Walter. I forgot to scroll over to see the whole plot.

[Image Analyst]

Did you simply try

f = fft(yourSignal);

And Emma, don't forget to click the link above to show the older comments.

Answer 1

[Meg Noah]

If I understand your question, it sounds like you want to find the underlying signal frequencies of your data.

clear variables
close all
clc
%% *bloodFlowRateSimulation*
% synthesized fake data of 3 frequencies is used to illustrate how to
% interpret fft results in matlab
%% *references*
% https://www.gaussianwaves.com/2015/11/interpreting-fft-results-complex-dft-frequency-bins-and-fftshift/
% https://blogs.uoregon.edu/seis/wiki/unpacking-the-matlab-fft/

rndStream = RandStream.create('mrg32k3a','Seed',6626231);

Fs = 1000;     % [samples/s] sampling frequency                    
T = 1/Fs;      % [s] sampling period       
N = 500;       % [samples] Length of signal
t = (0:N-1)*T; % [s] Time vector
deltaF = Fs/N; % [1/s]) frequency intervalue of discrete signal

% fake heart blood flow data is synthesis of 3 frequencies with 3 different
% amplitudes and phases
freqs  = [226 26 10];
amps   = [0.05 0.3 1];
dphase = [0 pi 0];
S = zeros(1,N);
for ifreq = 1:length(freqs)
    S = S + amps(ifreq)*sin(2*pi*freqs(ifreq)*t+dphase(ifreq));
end
X = S + abs(min(S(:))) + 0.1*(rand(rndStream,1,N)-0.5); 
 
figure('color','white','position',[70 100 600 900]); 
subplot(3,1,1);
plot(1e3*t(1:50),X(1:50))
title({'Signal + Zero-Mean Random Noise'; ...
    ['Signal Frequencies: ' num2str(freqs) ' Hz']})
xlabel('t (milliseconds)')
ylabel('X(t)')

% compute the fast fourier transform
Y = fft(X);

% manually shifting the FFT
% these values work for either odd or even sized arrays
Amp = [Y(ceil(end/2)+1:end) Y(1) Y(2:ceil(end/2))];
Amp = abs(Amp/N);
if (mod(N,2) == 0)
sampleIndex = -N/2:1:N/2-1; %raw index for FFT plot
else
sampleIndex = -(N-1)/2:1:(N-1)/2; %raw index for FFT plot
end

subplot(3,1,2);
plot(deltaF*sampleIndex, Amp);
hold on;
idx = find(Amp > 0.02);
plot(deltaF*sampleIndex(idx), Amp(idx), '+');
for k = 1:length(idx)
    if (idx(k) > (N-1)/2)
        text(deltaF*sampleIndex(idx(k)), Amp(idx(k))+0.15,...
            ['f=' num2str(deltaF*sampleIndex(idx(k))) ' Hz']);
    end
end
xlabel('Frequency [Hz]');
ylabel('Amplitude');
title('Interpeting data indices');

subplot(3,1,3);
half_f = deltaF*(0:(N/2));
plot(fftshift([half_f -fliplr(half_f(2:end+mod(N,2)-1))]), ...
    abs(fftshift(Y)/N));
xlabel('Frequency [Hz]');
ylabel('Amplitude');
title('Using fftshift');

This plot shows how to interpret the results. If there is a signal of a given frequency, you're accuracy in extracting it does depend on the time descritization. Is this the question you were asking, or did I completely misunderstand you?