Effacer les filtres
Effacer les filtres

length(dat​a.(Channel​First)) this worked for years, now it won't work WHY?

6 vues (au cours des 30 derniers jours)
Don
Don le 4 Mai 2022
Réponse apportée : Voss le 4 Mai 2022
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');

Réponses (2)

Steven Lord
Steven Lord le 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)

Voss
Voss le 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.

Catégories

En savoir plus sur EEG/MEG/ECoG dans Help Center et File Exchange

Produits


Version

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by