Biphasic Pulse signal with FFT

I want to create a biphasic pulse signal as used in Functional electrical stimulation, afterwards I would like to do a FFT to get a better understanding of which frequencies would be good for the stimulation. However my code has some errors. I changed my Input variables and now it says:

"Error using plot

Vectors must be the same length.

Error in BiphasicFFT (line 13)

plot(t,wholeSignal);"

What does it mean and how can I change that? I would like to change my parameters, so that I can analyze different stimulation parameters.

Is my FFT correct?

This is my code:

fs=100000; %sampling freq in Hz
pulsesPerS=60; %frequency in Hz
signalDuration=3; % in seconds
amplitude=22*10^-3; % in mA
pulseWidth=250*10^-6; % in µs

t=0:1/fs:signalDuration-1/fs; %vector over the duration of time and sampled every 1/fs

onePeriod=amplitude*[ones(round(pulseWidth*fs),1); -ones(round(pulseWidth*fs),1); zeros(round(1/pulsesPerS*fs)-2*round(pulseWidth*fs),1)];

wholeSignal=repmat(onePeriod,[pulsesPerS*signalDuration 1]);
figure (1)
plot(t,wholeSignal);

Error using plot
Vectors must be the same length.

xlabel('Time (seconds)');
ylabel('Amplitude');
title('Biphasic Pulse Train')

fourier=abs(fft(wholeSignal));
figure (2)
plot(fourier);
xlabel('Samples');
ylabel('Magnitude');
title('Frequency-Domain plot - FFT');

N=length(wholeSignal); %length of signal 
Freq=((0:N-1)/N)*fs; %frequency axis for FFT
Amp=abs(wholeSignal)/(N/2); %amplitude axis for FFT
figure (3)
plot(Freq,fourier);
%stem(Freq,fourier);
xlabel('Frequency in Hz');
ylabel('Amplitude');
title('FFT')

Answer 1

[Star Strider]

The problem was that ‘wholeSignal’ had 60 more elements than ‘t’ so I artificially shortened it:

wholeSignal = wholeSignal(1:numel(t));

so that they both fit.

After that, your code works. I would suggest doing a one-sided fft rather than two-sided, since it is a bit easier (more intuitive) to interpret. See the fft documentation for help with that.

Windowing the fft may also help, and you can do that with:

fourier=abs(fft(wholeSignal.*hann(numel(wholeSignal))));

That produces a slightly cleaner result. I added it (commented-out) so you can experiment with it.

fs=100000; %sampling freq in Hz
pulsesPerS=60; %frequency in Hz
signalDuration=3; % in seconds
amplitude=22*10^-3; % in mA
pulseWidth=250*10^-6; % in µs

t=0:1/fs:signalDuration-1/fs; %vector over the duration of time and sampled every 1/fs

onePeriod=amplitude*[ones(round(pulseWidth*fs),1); -ones(round(pulseWidth*fs),1); zeros(round(1/pulsesPerS*fs)-2*round(pulseWidth*fs),1)];

wholeSignal=repmat(onePeriod,[pulsesPerS*signalDuration 1]);
wholeSignal = wholeSignal(1:numel(t));

figure (1)
plot(t,wholeSignal);
xlabel('Time (seconds)');
ylabel('Amplitude');
title('Biphasic Pulse Train')

fourier=abs(fft(wholeSignal));
% fourier=abs(fft(wholeSignal.*hann(numel(wholeSignal))));
figure (2)
plot(fourier);
xlabel('Samples');
ylabel('Magnitude');
title('Frequency-Domain plot - FFT');

N=length(wholeSignal); %length of signal 
Freq=((0:N-1)/N)*fs; %frequency axis for FFT
Amp=abs(wholeSignal)/(N/2); %amplitude axis for FFT
figure (3)
plot(Freq,fourier);
%stem(Freq,fourier);
xlabel('Frequency in Hz');
ylabel('Amplitude');
title('FFT')

Those were the only changes I made to your code.

Good luck with your FES research!

.