load Parkinson_Demo_Data.mat

MEG=data.trial{1}(1,:);
EMG=data.trial{1}(2,:);
sf=data.fsample;

tmpTime=data.time;
tmpData=EMG;
Nsamples=size(tmpData,2);% if all dimentions but one are 1 then Nsamples=length(EMG,2);



%% EXERCISE 1
edit fftdemo %to execute commands from within fftdemo select them and press F9


%% EXERCISE 2
% QUESTION 1
Nfft=pow2(nextpow2(Nsamples));
Y = fft(tmpData,Nfft);
Pyy = Y.*conj(Y)/(Nfft);
f = sf*(0:(Nfft/2-1))/Nfft;
figure
plot(f,Pyy(1:Nfft/2));
title('Power spectral density')
xlabel('Frequency (Hz)')
% Things are not very clear because low frequencies seem to have much
% larger values than higher frequencies. When such big differences occur
% try plotting in db by taking 10*log10(X)
figure
logPyy=10*log10(Pyy)
plot(f,logPyy(1:Nfft/2));
title('Power spectral density(dB)')
xlabel('Frequency (Hz)')



% QUESTION 2 
estimator1=spectrum.periodogram;
Nfft=pow2(nextpow2(Nsamples));
figure
psd(estimator1,tmpData,'Fs',sf,'NFFT',Nfft)
% Check the effect of increasing NFFT
for i=8:nextpow2(Nsamples)
    i
    Nfft=pow2(i);
    psd(estimator1,tmpData,'Fs',sf,'NFFT',Nfft);
    pause
end



% QUESTION 3
Nfft=pow2(nextpow2(Nsamples));
estimator2=spectrum.welch;
estimator2 % OBSERVE THE PARAMETERS OF YOUR ESTIMATOR
set(estimator2,'WindowName') % CHECK THE LIST OF WINDOW TYPES AVAILABLE

tmpWinLength=floor(Nfft/10); % DO NOT FORGET THAT FFT USES NFFT DATA POINTS AND NOT THE ORIGINAL Nsamples OF YOUR DATA
set(estimator2,'SegmentLength',tmpWinLength)
figure
psd(estimator2,tmpData,'Fs',sf,'NFFT',Nfft);
% Check the effect of reducing the Window Length and consequently
% increasing the number of segments used for deriving the averaged
% spectrum.
for i=2:4:22
    tmpWinLength=floor(Nfft/i); % DO NOT FORGET THAT FFT USES NFFT DATA POINTS AND NOT THE ORIGINAL Nsamples OF YOUR DATA
    set(estimator2,'SegmentLength',tmpWinLength)
    psd(estimator2,tmpData,'Fs',sf,'NFFT',Nfft);
    pause
end
% The smaller the window and consequently the more the number of segments
% the data is split for averaging the spectrum, the lower the variance.


% Check the efect of using different types of Windows
Nfft=pow2(14);
tmpWinLength=floor(Nfft/3);
set(estimator2,'SegmentLength',tmpWinLength)
figure
for tmpWinName={'Hamming','Gaussian','Chebyshev','Flat Top'},
   tmpWinName
   set(estimator2,'WindowName',tmpWinName{:})
   psd(estimator2,tmpData,'Fs',sf,'NFFT',Nfft);
   pause
end
% Variance appears to become smaller as progressing in the above list of
% Window types. Check with wintool the amplitude characteristics of these
% window types. The spread of the main lobe seems to correlate with
% behavior of variance. Variance reduction is not very significant though.
% Typically a Hamming Window is the default choice for most applications.




% If you wanna have a go with a MultiTaper Method try the following:
% (Multitapers are used mainly to avoid Bias that might occur when
% windowing (welch) methods as above are used due to the fact that some
% segments might create high erroneous spectrum components that will
% prevail in the spectrum average).
estimator3=spectrum.mtm(15);
Nfft=pow2(11);
figure
subplot(2,1,1)
psd(estimator1,tmpData,'Fs',sf,'NFFT',Nfft);
subplot(2,1,2)
psd(estimator3,tmpData,'Fs',sf,'NFFT',Nfft);



%% EXERCISE 3
Nfft=pow2(11);
tmpWinLength=floor(Nfft/3);
tmpOverLap=floor(tmpWinLength/2);
tmpX=EMG;
tmpY=MEG;
mscohere(tmpX,tmpY,tmpWinLength,tmpOverLap,Nfft,sf)


figure
for i=10:15
    i
    tmpWinLength=512;
    tmpOverLap=256;
    Nfft=pow2(i);
    mscohere(tmpX,tmpY,tmpWinLength,tmpOverLap,Nfft,sf)
    pause
end


