FFT of EEG data

This code is for doing FFT on eeg data to findout the spectral amplitudes of Fundamental frequency and its harmonics,H2,h3 and H4.

Fundamental frequency is 100Hz

time window is from 0milliseconds to 156 milliseconds.

Sampling frequency is 20000Hz

Wheni run the code the following error is shown.

Unable to perform assignment because the left and right sides have a different number of elements.

Error in ff (line 70)

sample(1:length(env)) = env'.*sample(1:length(env));

Following is the code

% The input file should be in a txt format and all the headers should be
% removed. The file should only contain the waveform and no other details
% and should be placed in aseparate folder in the matlab directory.     %%
clear all;clc
epoch = 156; %the total epoch duration in milliseconds
prestim = 0;   % pre-stimulus time ignoring the negative sign
FFTstart = 0; %starting time range to run FFT
FFTend = 130 ;  %ending time range to run FFT
fund = 100;
bin = 10;
%% Do not edit beyond this point
post = epoch - prestim; %post stimulus time
lowf0=fund-bin;  %F0 starting  frequncy
highf0=fund+bin; %F0 ending frequency
lowh2 = 2*fund - bin; %Second harmonic starting frequency
highh2 = 2*fund + bin; % Second harmonic ending frequency 
lowh3 = 3*fund - bin; %H3 starting frequency
highh3 = 3*fund + bin; % H3 ending frequency
lowh4 = 4*fund - bin; %H4 starting frequency
highh4 = 4*fund + bin; %H4 ending frequency
[path] = uigetdir('C:\Users\Gowri\Documents\MATLAB\*.txt', 'select the folder with the ASCII files');
path = [path '\*.txt'];
files = dir(path);
nfiles = length(files);
randlist=1:1:nfiles;
files = char(files.name);
a = length(files);
FFTstart = FFTstart/1000;
FFTend = FFTend/1000;
%epoch details
for i= 1:nfiles
    %number of points  and timescale
    name = files(i,:);
    a = dlmread(name);
    n = length(a);
    fs = n*1000/epoch;
    time = -prestim/1000 :1/fs: post/1000 -1/fs;
    zerotime = -prestim/1000:1/fs:0-1/fs ;
    %baseline correction ref: prestimulus
    detrend = a - mean(a(1:size(zerotime)));
    %%spectral analysis%%
    %Time slice%
    [trash position1] = min(abs(time-FFTstart));
    [trash position2] = min(abs(time-FFTend));
    sample = (a(position1:position2))';
    %zero padding%
    lin = fs;
    win = round(fs*(5/1000));
    window = hann(win);
    env = window(1:round((length(window))/2));
    vne = fliplr(env');
    sample(1:length(env)) = env'.*sample(1:length(env));
    sample(length(sample)-length(vne)+1:length(sample)) = vne.*sample(length(sample)-length(vne)+1:length(sample));
    zero = zeros(1,round(fs-length(sample)));
    abc = [sample zero];
    %FFT
    spect = abs(fft(abc));
    amp=spect.*(2./length(sample));
    amp = amp(1:round(round(fs)/2));
    n=length(spect);
    freq=fs/n.*(1:n);
    f=freq(1:(round(n/2)));
    %%H1
    fstartL =  find(f<=lowf0);
    fstartL = fstartL(1,length(fstartL));
    fstart = f(1,fstartL);
    fendL = find(f>=highf0);
    fendL = fendL(1,1);
    fend = f(1,fendL);
    range = f(fstartL:fendL);
    f0 = amp(fstartL:fendL);
    %find frequency with maximum amplitude
    MAX = max(f0);
    maxlocus = find(f0==MAX);
    Fo = range(1,maxlocus);
    H1 = mean(f0);
    % H2
    fstartL =  find(f<=lowh2);
    fstartL = fstartL(1,length(fstartL));
    fstart = f(1,fstartL);
    fendL = find(f>=highh2);
    fendL = fendL(1,1);
    fend = f(1,fendL);
    range = f(fstartL:fendL);
    h2 = amp(fstartL:fendL);
    H2 = mean(h2);
    %H3
    fstartL =  find(f<=lowh3);
    fstartL = fstartL(1,length(fstartL));
    fstart = f(1,fstartL);
    fendL = find(f>=highh3);
    fendL = fendL(1,1);
    fend = f(1,fendL);
    range = f(fstartL:fendL);
    h3 = amp(fstartL:fendL);
    H3 = mean(h3);
    %H4
    fstartL =  find(f<=lowh4);
    fstartL = fstartL(1,length(fstartL));
    fstart = f(1,fstartL);
    fendL = find(f>=highh4);
    fendL = fendL(1,1);
    fend = f(1,fendL);
    range = f(fstartL:fendL);
    h4 = amp(fstartL:fendL);
    H4 = mean(h4);
    Head = {'subject' 'F0' 'F0Amp' 'H2amp' 'H3amp' 'H4amp'};
    M = {name Fo H1 H2 H3 H4 };
    header = Head;
    analysis(i+1,:) = M;
end
analysis(1,:)=header;
%xlswrite('baseline.txt', analysis)

[Rik]

Without the data files it is not possible to reproduce your error, but I'm guessing env is a vector in a different direction than you expected, leading to this:

[1 2].*[3;4]

ans = 2×2

3 6 4 8

And obviously that will not fit in two elements.

You should be aware that the apostrophe is the complex conjugate, not just a transpose. The difference is non-existent for real data, but since you're working with fft, complex data might show up. You need .' instead. You should also reconsider your use of length. It doesn't do what you might want it to do. Did you actually intend to have max(size(___))? You probably want to use numel instead.

[Gowri R]

This is the eeg file which i used.

Iam a beginner in matlab and this code was written by somebody else. I was trying to analyze my file with this.

Comment 1: "You need .' instead " - this i tried once.same error image is showing up.

comment 2: "You should also reconsider your use of length. It doesn't do what you might want it to do. Did you actually intend to have max(size(___))? You probably want to use numel instead" - this i couldnt understand.

Could you please simplify these things for me if possible

Answer 1

[Rik]

Are you sure your file is in the format this function expects?

This is a good lesson for your own code: make sure to write documentation for your functions that explain what kind of data is required and in what format. There are hardly any comments in this code, making debugging very difficult.

I didn't mean that replacing the complex conjugate with a transpose would solve your issue, I merely wanted to alert you to the difference between the two, which may be very important in some contexts.

Regarding your last question: the length function should almost never be used. Calling length(A) is equivalent to calling max(size(A)), which is a very odd thing to want. Generally when people use the length function they actually want the total number of elements in an array, which is what the numel function will provide.

If you have trouble with Matlab basics you may consider doing the Onramp tutorial (which is provided for free by Mathworks). Make sure to read the documentation for any function you do not understand or which is giving you results you don't understand.

But now on to your main question. Let's see how far this code gets. (we need to edit the part where it selects a local folder on your computer)

% The input file should be in a txt format and all the headers should be
% removed. The file should only contain the waveform and no other details
% and should be placed in aseparate folder in the matlab directory.     %%
clear all;clc
epoch = 156; %the total epoch duration in milliseconds
prestim = 0;   % pre-stimulus time ignoring the negative sign
FFTstart = 0; %starting time range to run FFT
FFTend = 130 ;  %ending time range to run FFT
fund = 100;
bin = 10;
%% Do not edit beyond this point
post = epoch - prestim; %post stimulus time
lowf0=fund-bin;  %F0 starting  frequncy
highf0=fund+bin; %F0 ending frequency
lowh2 = 2*fund - bin; %Second harmonic starting frequency
highh2 = 2*fund + bin; % Second harmonic ending frequency 
lowh3 = 3*fund - bin; %H3 starting frequency
highh3 = 3*fund + bin; % H3 ending frequency
lowh4 = 4*fund - bin; %H4 starting frequency
highh4 = 4*fund + bin; %H4 ending frequency
%[path] = uigetdir('C:\Users\Gowri\Documents\MATLAB\*.txt', 'select the folder with the ASCII files');
%path = [path '\*.txt'];
files = dir('*.txt');
nfiles = length(files);
randlist=1:1:nfiles;
files = char(files.name);
a = length(files);
FFTstart = FFTstart/1000;
FFTend = FFTend/1000;
%epoch details
for i= 1:nfiles
    %number of points  and timescale
    name = files(i,:);
    a = dlmread(name);
    n = length(a);
    fs = n*1000/epoch;
    time = -prestim/1000 :1/fs: post/1000 -1/fs;
    zerotime = -prestim/1000:1/fs:0-1/fs ;
    %baseline correction ref: prestimulus
    detrend = a - mean(a(1:size(zerotime)));
    %%spectral analysis%%
    %Time slice%
    [trash position1] = min(abs(time-FFTstart));
    [trash position2] = min(abs(time-FFTend));
    sample = (a(position1:position2))';
    %zero padding%
    lin = fs;
    win = round(fs*(5/1000));
    window = hann(win);
    env = window(1:round((length(window))/2));
    vne = fliplr(env');
    
    fprintf('native sizes of env and sample:\n')
    size(env),size(sample)
    fprintf('sizes of env.'' and sample(1:numel(env)):\n')
    size(env'),size(sample(1:length(env)))
    sample(1:length(env)) = env'.*sample(1:length(env));
    % skip the remainder for brevity of the post
end

native sizes of env and sample:

ans = 1×2

50 1

ans = 1×2

2605 1

sizes of env.' and sample(1:numel(env)):

ans = 1×2

1 50

ans = 1×2

50 1

Unable to perform assignment because the left and right sides have a different number of elements.

As you can see, you are multiplying a vector of 1x50 with a vector of 50x1, which will implicitly expand to 50x50. Removing the conjugation will solve this problem.

The original author should have put in some checks to confirm the dimensions of the input arrays.