length(data.(ChannelFirst)) this worked for years, now it won't work WHY?
    5 Ansichten (letzte 30 Tage)
  
       Ältere Kommentare anzeigen
    
Select edf data file to analyze
Enter SAMPLE RATE for this subject(check log book!):    512
Enter name of first channel to analyze :  F5
Enter name of second channel to analyze :  F3
Unrecognized field name "F5".
Error in SpectralRatioST1 (line 162)
veclength = length(data.(ChannelFirst));         % length of input EEG vector
Here is the input file:
tried to put it in but no luck
its an EDF file
Here is the program:
% SpectralRatioST -- read multi-channel file of EEG traces, calculate power
%        spectra,form ratio of frequency bands (e.g., Theta/low Alpha), and plot
%        ratio over time
% This version uses Structured Array
%
% Author: Don Scott, 2012.  
clear;
clf;
clc;
% Read EEG file data
% % ask user
display('Select edf data file to analyze');
[filename, filepath] = uigetfile('*.edf', 'Select edf file'); 
if filename == 0 return; end;
fname = [filepath filename];
% ask user
cd        'C:\Projects\ParkinsonsProject\Paul Favko\'
[hdr, record]=edfRead(fname);
if length(hdr.label)< 50 
    data.C4 = record(6,:);
    data.P3 = record(7,:);
    data.P4 = record(8,:);   
    data.O1 = record(9,:);
    data.O2 = record(10,:);
    data.F7 = record(11,:);
    data.F8 = record(12,:);
    data.T3 = record(13,:);
    data.T4 = record(14,:);
    data.T5 = record(15,:);
    data.T6 = record(16,:);
    data.A1 = record(17,:);
    data.A2 = record(18,:);
    data.Fz = record(19,:);
    data.Cz = record(20,:);
    data.Pz = record(21,:);
    data.SubL = record(22,:);
    data.SubR = record(23,:);
    data.ECG = record(24,:);
else
    %        %data.Fp1 = record(1,:);
    %        %data.Fp2 = record(2,:);
    %        data.F4 = record(3,:);
    %        data.C3 = record(4,:);
    %        data.C4 = record(5,:);
    %        data.P3 = record(6,:);
    %        data.P4 = record(7,:);
    %        %data.O1 = record(8,:);
    %        %data.O2 = record(9,:);
    %        %data.F7 = record(10,:); 
    %        %data.F8 = record(11,:);
    %        data.FC3 = record(12,:);
    %        %data.FT7 = record(13,:);
    %        %data.FT8 = record(14,:);
    %        data.T7 = record(15,:);
    %        data.T8 = record(16,:);
    %        data.A1 = record(17,:);
    %        data.A2 = record(18,:);
    %        data.FCz = record(19,:);
    %        data.TP7 = record(20,:);
    %        data.CPz = record(21,:);
    %        data.CP3 = record(22,:);
    %        %data.P7 = record(23,:);
    %        data.TP8 = record(24,:);
    %        %data.P8 = record(25,:);
    %        data.CP4 = record(26,:);
    %        %data.Oz = record(27,:);
    %        %data.HEOL = record(28,:);
    %        %data.HEOR = record(29,:);
    %        %data.FPz = record(30,:); 
    %        %data.AF3 = record(31,:);
    %        %data.AF7 = record(32,:);
    %        data.F5 = record(33,:);
    %        %data.AF8 = record(34,:);
    %        %data.AF4 = record(35,:);
    %        data.F1 = record(36,:);
    %        data.FC5 = record(37,:);
    %        data.F6 = record(38,:);
    %        data.F2 = record(39,:);
    %        data.FC1 = record(40,:);   
    %        data.C5 = record(41,:);
    %        data.FC6 = record(42,:);
    %        data.FC2 = record(43,:);
    %        data.C2 = record(44,:);
    %        data.C1 = record(45,:);
    %        data.CP1 = record(46,:);
    %        data.CP5 = record(47,:);
    %        %data.P5 = record(48,:);
    %        %data.PO7 = record(49,:);
    %        %data.PO8 = record(50,:); 
    %        data.C6 = record(51,:);
    %        data.CP6 = record(52,:);
    %        %data.PO6 = record(53,:);
    %        %data.P6 = record(54,:);
    %        data.CP2 = record(55,:);
    %        %data.PO4 = record(56,:);
    %        %data.P2 = record(57,:);
    %        %data.PO2 = record(58,:);
    %        %data.P1 = record(59,:);
    %        %data.PO3 = record(60,:); 
    %        %data.PO5 = record(61,:); 
    %        data.Fz=record(39,:);    %duplicated because cap shows 'Fz' but header
    %                                 % shows 'F2'. rest of program uses 'Fz'
    %                                 % instead of T3,T4  use T7,T8
    %        data.F3=record(40,:);
    %        data.Cz=record(41,:);
    %        data.FC4=record(42,:);
    data.F3 = record(1,:);
    data.F4 = record(2,:);
    data.C3 = record(3,:);
    data.C4 = record(4,:);
    data.P3 = record(5,:);
    data.P4 = record(6,:);
    data.CPz = record(7,:);
    data.CP3 = record(8,:);
    data.Fz = record(9,:);
    data.Cz = record(10,:);
    data.CP4 = record(11,:);
    data.FC3 = record(12,:);
    data.TP8 = record(13,:);
    data.TP7 = record(14,:);
    data.FCz = record(15,:);
    data.FC4 = record(16,:);
    data.F5 = record(17,:);
    data.F1 = record(18,:);
    data.FC5 = record(19,:);
    data.F6 = record(20,:);
    data.F2 = record(21,:);
    data.FC1 = record(22,:);   
    data.C5 = record(23,:);
    data.FC6 = record(24,:);
    data.FC2 = record(25,:);
    data.C2 = record(26,:);
    data.C1 = record(27,:);
    data.CP1 = record(28,:);
    data.CP5 = record(29,:);
    data.C6 = record(30,:);
    data.CP6 = record(31,:);
    data.CP2 = record(32,:);
    % data.Fz = record(33,:);    %duplicated because cap shows 'Fz' but header
    % shows 'F2'. rest of program uses 'Fz'
    % instead of T3,T4  use T7,T8
    %data.F3 = record(34,:);
    %data.Cz = record(35,:);
    %data.FC4 = record(36,:);
end;
srate = input ('Enter SAMPLE RATE for this subject(check log book!):    ');   
ChannelFirst = input ('Enter name of first channel to analyze :  ','s');
ChannelSecond = input ('Enter name of second channel to analyze :  ','s');
veclength = length(data.(ChannelFirst));         % length of input EEG vector
iteration=0;
lb1secint = 0;                                   % lower bound of 1-sec segment of input file 
ub1secint = 0;                                   % upper bound of 1-sec segment of input file 
RatioResults=[0 0 0 0 0 0 0 0 0]';  % initialize matrix of ratios X sec for output plot
RatioResultsID = ['seconds  ' 'RatioTheta1LoAlpha1  ' 'RatioTheta1HiAlpha1  ' ...
    'RatioTheta2LoAlpha2  ' 'RatioTheta2HiAlpha2 ' 'Theta1 ' 'Theta2 ' ...
    'HiAlpha1 ' 'HiAlpha2 '];
HiSpike=[0 0 0];  % initialize matrix of ratios X sec for Theta/HiAlpha Spikes
HiSpikeID = ['seconds  ' 'Channel Fz  ' 'Channel P3  '];     %'Theta1  '  'HiAlpha1  ' ...
%  'Theta2'
%  'HiAlpha2']
% sample rate = 512, take input vector in 1-sec increments
ub1secint = (iteration+1)*srate;    % lower bound of 1-sec segment of input file 
lb1secint = (iteration*srate)+1;    % upper bound of 1-sec segment of input file 
disp(' Diagnostic values will report spectral details of Theta, HiAlpha ');
disp(' and Theta/HiAlpha ratio for selected channels.  Set up interation ');
disp(' limits in lines 130 and 131 of SpectralRatio.m');
DiagnosticChoice = input(' Do you want  diagnostic values for this file? (y/n) ','s');
j=1;                        %initialize index
iteration=1;
while ub1secint <= veclength;   
    fprintf ('iteration  %6.0f\n',iteration);                                
    %fprintf ('LINE 75 low bound(lb1secint): %6.0f   up bound(ub1secint): %6.0f\n', lb1secint,ub1secint);
    fnme1=data.(ChannelFirst)(lb1secint:1:ub1secint);  % vecvtor of current 1-sec interval, first channel selected
    fnme2=data.(ChannelSecond)(lb1secint:1:ub1secint);  % vecvtor of current 1-sec interval, second channel selected
    %fprintf ('length of fnme1: %6.0f\n', length(fnme1));
    %fprintf ('length of fnme2: %6.0f\n', length(fnme2));
    %fprintf('  \n\n');
    % Calculate Power Spectra
    %x=spect1.f(535:1:750);     %get freq range 0 to 50 Hz
    %y=spect1.P(535:1:750);     %get uv^2 in range 0 to 50 Hz
    [spectra1, freq] = pwelch(fnme1,hamming(length(fnme1)),[],2*length(fnme1),srate);
    power1=10*log10(spectra1);
    [spectra2, freq] = pwelch(fnme2,hamming(length(fnme2)),[],2*length(fnme2),srate);
    power2=10*log10(spectra2);
    % Form Frequency Band ratio  Theta 4.5-8 Hz, HiALPHA  11.5 - 14 Hz
    Theta1   =mean(power1(9:16));
    LoAlpha1 =mean(power1(17:20));
    HiAlpha1 =mean(power1(23:28));
    Theta2   =mean(power2(9:16));
    LoAlpha2 =mean(power2(17:20));
    HiAlpha2 =mean(power2(23:28));
    RatioTheta1LoAlpha1 = Theta1/LoAlpha1;
    RatioTheta1HiAlpha1 = Theta1/HiAlpha1;
    RatioTheta2LoAlpha2 = Theta2/LoAlpha2;
    RatioTheta2HiAlpha2 = Theta2/HiAlpha2;
    %Collect date from this iteration for additiion to output array & file
    thispass = [iteration RatioTheta1LoAlpha1 RatioTheta1HiAlpha1 ...
        RatioTheta2LoAlpha2 RatioTheta2HiAlpha2 Theta1 Theta2 HiAlpha1 HiAlpha2]';
    RatioResults=[RatioResults thispass];
    %fprintf('thispass  %4.2f\n',thispass)
    %  Diagnostic code to display Theta,HiAlpha spectral data
    %  for selected iterations
    %if iteration ==126;
    %fprintf('RatioTheta1HiAlpha1  %6.3f\n',RatioTheta1LoAlpha1);
    %fprintf('RatioTheta2HiAlpha2  %6.3f\n',RatioTheta2LoAlpha2);
    %break;
    %end;
    switch DiagnosticChoice
        case 'y'
            %if RatioTheta2HiAlpha2=='NaN';
            %fprintf('RatioTheta1HiAlpha1  %6.3f\n',RatioTheta1LoAlpha1);
            %end;
            %if RatioTheta1HiAlpha1=='NaN';
            %fprintf('RatioTheta2HiAlpha2  %6.3f\n',RatioTheta2LoAlpha2);
            %end;
            %if iteration ==127;
            %   break;
            %end;
            if iteration < 200
                precision=3;
                %  figure (6);
                % hold on;
                vecc=num2str(power1(9:16)', precision);
                Ch1out=['SpecThetaVals ' vecc];
                Theta1out=['Theta1 ',num2str(Theta1)];
                vecc=num2str(power1(23:28)', precision);
                Ch1out=['SpecHAlphaVals ' vecc];
                % plot(power1(23:28));
                HiAlpha1out=['HiAlpha1 ',num2str(HiAlpha1)];
                RatioTheta1HiAlpha1out = ['RatioTheta1HiAlpha1 ',num2str(RatioTheta1HiAlpha1)];
                if abs(RatioTheta1HiAlpha1) > 20 
                    %if iteration ==126
                    disp('1stChannel:  ');  %info detail for Channel 1 spectra, ratio
                    disp([Ch1out]); 
                    disp([Theta1out]);
                    disp([Ch1out]);
                    disp([HiAlpha1out]);
                    disp([RatioTheta1HiAlpha1out]);
                    disp('-----------------------------------------------------------');
                end;
                vecc=num2str(power2(9:16)', precision);
                Ch2out=['SpecThetaVals ' vecc];
                Theta2out=['Theta2 ',num2str(Theta2)];
                vecc=num2str(power2(23:28)', precision);
                Ch2out=['SpecHA2phaVals ' vecc];
                HiAlpha2out=['HiAlpha2 ',num2str(HiAlpha2)];
                RatioTheta2HiAlpha2out = ['RatioTheta2HiAlpha2 ',num2str(RatioTheta2HiAlpha2)];
                if abs(RatioTheta2HiAlpha2) > 20
                    disp('2ndChannel:  ');          %info detail for Channel 2 spectra, ratio
                    disp([Ch2out]); 
                    disp([Theta2out]);
                    disp([Ch2out]); 
                    disp([HiAlpha2out]);
                    disp([RatioTheta2HiAlpha2out]);
                    disp('===============================================================');
                end                                     %  hold off;
                %         end;
            end;
    end;
    % load Theta HiAlpha values into array HiSpike
    %HiSpike(1,iteration) = iteration;
    %HiSpike(4,iteration) = Theta1;    %value of mean Theta1
    %HiSpike(5,iteration) = HiAlpha1;  %value of mean HiAlpha1
    %HiSpike(6,iteration) = Theta2;    %value of mean Theta2
    %HiSpike(7,iteration) = HiAlpha2;  %value of mean HiAlpha2
    %  set boundaries for next 1-sec interval
    ub1secint = (iteration+1)*srate;
    lb1secint = (iteration*srate)+1;
    iteration=iteration+1;  
    j=j+1;
end
disp('after loop')
% Plot Ratio, write output files
ssRR3=isnan(RatioResults(3,:));
ssRR5=isnan(RatioResults(5,:));
RR3=RatioResults(3,:);
RR5=RatioResults(5,:);
stdevT1HiA1 = std(RR3(~ssRR3));
stdevT2HiA2 = std(RR5(~ssRR5));
%fprintf('stdevs %6.3f  %6.3f\n',stdevT1HiA1, stdevT2HiA2);
plot(RatioResults(1,:),RatioResults(3,:),'r-','LineWidth',2);
hold on
plot(RatioResults(1,:),RatioResults(5,:),'b-','LineWidth',2);
% Plot Events and noise
Event1      = [2 2;-20 30];    %Locate Button Press Event
Event2      = [152 152;-20 30];  %Button Press
%Event3      = [153 153;-20 30];  %Button Press
%Event4      = [72 72;-20 30];  %Button Press
%Event5      = [98 98;-20 30];  %Button Press
%Event6      = [109 109;-20 30];  %Button Press
%Event7      = [129 129;-20 30];  %Button Press
%Event8      = [157 157;-20 30];  %Button Press
%Event9      = [170 170;-20 30];  %Button Press
%Event10     = [182 182;-20 30];  %Button Press
%plot(Event1(1,:), Event1(2,:));
%plot(Event2 (1,:), Event2(2,:));
%plot(Event3 (1,:), Event3(2,:));
%plot(Event4 (1,:), Event4(2,:));
%plot(Event5 (1,:), Event5(2,:));
%plot(Event6 (1,:), Event6(2,:));
%plot(Event7 (1,:), Event7(2,:));
%plot(Event8 (1,:), Event8(2,:));
%plot(Event9 (1,:), Event7(2,:));
%plot(Event10 (1,:), Event8(2,:));
%plot(Event0 (1,:), Event3(2,:),'g--', 'LineWidth',2);
%Noise1 = [5 5; -30 30]; %locate beginning of Noise
%Noise1a = [21 21; -30 30]; %locate end of Noise
%Noise2 = [56 56; -20 30]; %locate beginning of Noise
%Noise2a = [109 109; -20 30]; %locate end of Noise
%Noise3 = [142 142; -20 30]; %locate beginning of Noise
%Noise3a = [143 143; -20 30]; %locate end of Noise
%Noise4 = [119 119; -20 30]; %locate beginning of Noise
%Noise4a = [120 120; -20 30]; %locate end of Noise
%plot(Noise1(1,:), Noise1(2,:),'r--','LineWidth',2);
%plot(Noise1a(1,:), Noise1a(2,:),'r--','LineWidth',2);
%plot(Noise2(1,:), Noise2(2,:),'r--','LineWidth',2);
%plot(Noise2a(1,:), Noise2a(2,:),'r--','LineWidth',2);
%plot(Noise3(1,:), Noise3(2,:),'r--','LineWidth',2);
%plot(Noise3a(1,:), Noise3a(2,:),'r--','LineWidth',2);
%plot(Noise4(1,:), Noise4(2,:),'r--','LineWidth',2);
%plot(Noise4a(1,:), Noise4a(2,:),'r--','LineWidth',2);
hold off
mytitle=['Plot of  ' filename, ' ', ChannelFirst,' ', ChannelSecond];
title(mytitle); 
%legend('Fz','P3','MWEvents','EEGNoise','Location','SouthWest');
%legend('Fz','P3','MWEvents','Location','SouthWest');
%legend(ChannelFirst, ChannelSecond,'MW Events','Location','SouthWest');
legend(ChannelFirst, ChannelSecond,'Location','SouthWest');
axis ([0 200 -250 250]);
xlabel ('Seconds');
ylabel('Theta/HiAlpha');
%Saving Theta/Alpha data files RatioResults
[p,n,e]=fileparts(filename);
%newFileName2 = fullfile(filepath, [' Channels',' ',ChannelFirst,' ',ChannelSecond]);
%save(newFileName, 'RatioResults','-ascii');
RatioResults(:,1)=[];  %remove first column of 0's
newFileName = fullfile(pathname,[n,' RatioResults',' ',ChannelFirst,' ',ChannelSecond,'.xls']);
%try
%    xlswrite(newFileName2,ChannelFirst);
%    xlswrite(newFileName2,ChannelSecond);
%catch
%    pause(3);
%    xlswrite(newFileName2,ChannelFirst);
%    xlswrite(newFileName2,ChannelSecond);
%end;
% try
% xlswrite(newFileName,RatioResults);
% catch
%     pause(5);
%     xlswrite(newFileName,RatioResults);
% end
try
    writematrix(RatioResults,newFileName);
catch
    pause(5);
    writematrix(RatioResults,newFileName);
end
%This section tests Theta/HiAlpha ratio for magnitude > 2*std Dev, 
%if true writes value into array HiSpike. Does for ChannelFirst 
%and ChannelSecond
% Feb 27 2013 problem with extreme values of Theta/HiAlpha, occur because
% HIAlpha gets very close to zero.  Switch cutoff to depend on second 
% quartile (= median)
lengthT1HiA1 = length(RatioResults(3,:));
lengthT2HiA2 = length(RatioResults(5,:));
cutoffSD = 1.96*stdevT1HiA1;
%cutoffQ=quantile(RatioResults(3,:),.99);
cutoffQ=quantile(RR3(~ssRR3),.99);
fprintf('CutoffSD = %6.3f  CutoffQ99 = %6.3f.  (QUANTILE99 used here)\n', cutoffSD,cutoffQ);
for i=1:lengthT1HiA1
    HiSpike(i,1)=RatioResults(1,i);
    if abs(RatioResults(3,i)) > cutoffQ;
        HiSpike(i,2)=RatioResults(3,i);
    end
    if abs(RatioResults(5,i)) > cutoffQ;
        HiSpike(i,3)=RatioResults(5,i);
    end;
end
%Saving Theta/Alpha data files HiSpike
newFileName1 = fullfile(pathname, [n,' HiSpikeST',' ',ChannelFirst,' ',ChannelSecond]);
try
    xlswrite(newFileName1,HiSpike);
catch
    pause(10);
    xlswrite(newFileName1,HiSpike);
end
% use input stmt here to get names for channels in ratios
%  format is     xyz=input('string');
%concatenate to get output file name
%outputfile = ['ThetaHiAlpha' ' filename' ' fname1' ' fname2' '-xls')
fprintf ('     Files Saved\n');
0 Kommentare
Antworten (2)
  Steven Lord
    
      
 am 4 Mai 2022
        Set an error breakpoint and run your code. When MATLAB stops on that line, show us the output of the following command:
whos data ChannelFirst
As stated on this documentation page "The dynamic fieldname can return either a character vector or a string scalar." My guess is that your code used to return a character vector or a string scalar but now returns something else (a cell array containing character vectors, which can be called a cellstr: see the iscellstr function, or a non-scalar string) and the output of whos should tell us that.
If that indicates that ChannelFirst is a character vector or string scalar then can you show us what fields the data variable has (assuming it's a struct array?)
fieldnames(data)
0 Kommentare
  Voss
      
      
 am 4 Mai 2022
        The error says "F5" is not a field of the struct data.
Looking at the code, I can see that data.F5 doesn't get assigned if length(hdr.label)< 50.
If I were you, I would inspect the hdr and record returned from edfread, and make sure they are what you expect.
0 Kommentare
Siehe auch
Kategorien
				Mehr zu EEG/MEG/ECoG finden Sie in Help Center und File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!


