length(data.(ChannelFirst)) this worked for years, now it won't work WHY?
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Respuestas (2)
Steven Lord
el 4 de Mayo de 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 comentarios
Voss
el 4 de Mayo de 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 comentarios
Ver también
Categorías
Más información sobre EEG/MEG/ECoG en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!