Create a new version of paste:

function plot_Fourier1(inputSignal,fs,wlen)

% fs = sampling frequency
% inputSignal = single column vector of numerical data
% wlen = window length (recomended to be an integer power of 2)
% shorter windows give higher resolution and are less affected by
% high-frequency noise, they do however tend to loose the windowing effect
% if they are chosen too small wrt the length of inputSignal
 
% To run the function on the desired EMG data, run this from the command
% line:
% >> plot_Fourier1(Kanaal1mVolt,1000,512); (don't forget to import
% Kanaal1mVolt)


signal = inputSignal;

x = signal(1: (length(signal)-1)); % exclude last value (NaN) because signal processor can't deal with that...
time = 1 : (length(signal)-1);
time = time.';
time = time/1000; %divide by sampling frequency

%%%%%%%%%%%%%%%%%%%%%%Plot Raw Data%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%define analysis parameters for stft
xlen = length(x);    % length of the signal 
hop = wlen/4;       % hop size (recomended to be power of 2)
nfft = 2^nextpow2(xlen); 

figure(1);
plot(time,x);
set(gca, 'FontName', 'Arial', 'FontSize', 14);
xlabel('Time, s');
ylabel('Amplitude, mV');
title('Time Domain Signal - Raw EMG');

% apply short-time Fourier transform for spectral analysis

% number of fft points (recomended to be power of 2)
% perform STFT
[S, f, t] = stft(x, wlen, hop, nfft, fs);

K = sum(hamming(wlen, 'periodic'))/wlen;
% take the amplitude of fft(x) and scale it, so not to be a
% function of the length of the window and its coherent amplification
S = abs(S)/wlen/K;

if rem(nfft, 2)                     
    S(2:end, :) = S(2:end, :).*2;
else                                
    S(2:end-1, :) = S(2:end-1, :).*2;
end

% convert amplitude spectrum to dB (min = -120 dB)
S = 20*log10(S + 1e-6);

% plot the spectrogram

figure(2);
surf(t, f, S);
colormap(jet);
shading interp;
axis tight;
box on;
view(0, 90);
handl = colorbar;
set(handl, 'FontName', 'Arial', 'FontSize', 14,'Colormap',jet);
xlabel('Time, s');
ylabel('Frequency, Hz');
title('Spectogram of Raw EMG');
ylabel(handl, 'Magnitude, dB');


%%%%%%%%%%%%%%%%%%Filter the Data%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%notch filter to remove noise from 50 Hz
% Design a filter with a Q-factor of Q=35 to remove a 50 Hz tone from 
% system running at 1000 Hz.
Wo = 50/(fs/2);  BW = Wo/35;
[b,a] = iirnotch(Wo,BW);  
x = filter(b,a,x); %use filtfilt to avoid phase shifts?

% Get rid of the DC component (hp filter at 350 Hz)
fc=350;
[b,a] = butter(6,2*fc/fs,'high');
x = filter(b,a,x);


%filter out noise from the electrical setup, including higher harmonics
%of the AC supply (aka evrthg over 450 Hz)
fc=450;
[b,a] = butter(6,2*fc/fs,'low');
x = filter(b,a,x);


% number of fft points (recomended to be power of 2)
% perform STFT
[S1, f, t] = stft(x, wlen, hop, nfft, fs);
% define the coherent amplification of the window
%'periodic' â€” This option is useful for spectral analysis because it
% enables a windowed signal to have the perfect periodic extension implicit 
% in the discrete Fourier transform. When 'periodic' is specified, hamming 
% computes a window of length L + 1 and returns the first L points.

K = sum(hamming(wlen, 'periodic'))/wlen;
% take the amplitude of fft(x) and scale it, so not to be a
% function of the length of the window and its coherent amplification
S1 = abs(S1)/wlen/K;


if rem(nfft, 2)                   
    S1(2:end, :) = S1(2:end, :).*2;
else                               
    S1(2:end-1, :) = S1(2:end-1, :).*2;
end

% convert amplitude spectrum to dB (min = -120 dB)
S1 = 20*log10(S1 + 1e-6);

% plot the spectrogram
figure(3);
surf(t, f, S1);
colormap(jet);
shading interp;
axis tight;
box on;
view(0, 90);
handl = colorbar;
set(handl, 'FontName', 'Arial', 'FontSize', 14, 'Colormap', jet);
xlabel('Time, s');
ylabel('Frequency, Hz');
title('Spectrogram of Filtered EMG');
ylabel(handl, 'Magnitude, dB');
set(gca,'Ylim',[0 550]);