Why is my plot showing up blank?

I am attempting to plot a frequency spectrum for some MEG data using this script seedssE.m, originally written by Nai Ding for his 2016 study looking at cortical tracking of linguistic structures: https://www.nature.com/articles/nn.4186

clear
load('R1167_test.mat')
for trl=1:size(data,3)
    x(:,:,trl)=resample(data(:,:,trl),1,1);end
x=single(x);
save R1167-1000Hz_test x
 
%%
clear;
load('R1167-1000Hz_test.mat')
fs=1000;
[b]=fir1(200,[0.5 6]/(fs/2));a=1;
clear xf
for trl=1:size(x,3)
    xf(:,:,trl)=filter(b,a,x(:,:,trl));
end
xf=xf(fs+[1:0.32*fs*4*10],:,:);
 
x0=unfold(xf);
x1=mean(xf,3);
c0=x0'*x0;
c1=x1'*x1;
 
keep2=10.^-12;keep1=[];
[todss,fromdss,ratio,pwr]=dss0(c0,c1,keep1,keep2);
 
y=fold(unfold(xf)*todss(:,1:10)*fromdss(1:10,:),size(xf,1));
[sa,sb,f]=plot_itc_avg_spectrum(y,fs,[0 5]);
subplot 211;set(gca,'xtick',[1/4 1/2 1 2 4]/0.32);grid
hold on;plot(f,mean(sa,2)*3,'k','linewidth',2)
subplot 212;set(gca,'xtick',[1/4 1/2 1 2 4]/0.32);grid
hold on;plot(f,mean(sb,2)*3,'k','linewidth',2)

This makes use of the function plot_itc_avg_spectrum, shown in full below:

function [itc_sprectrum,avg_spectrum,f]=plot_itc_avg_spectrum(eeg_data,fs,display_frequency_range,f_label)
%     calculate and plot the neural response spectrum
%     input:
%         eeg_data (time * channel * trial)
%         fs (sampling rate in Hz)
%         display_frequency_range (the frequency range that will be plotted)
%         f_label (optional, frequencies being labeled)
%     output:
%         itc_sprectrum (frequency * channel, inter-trial phase coherence spectrum R^2),
%         avg_spectrum (frequency * channel, power spectrum of the response averaged over trials)
%         f (frequency labels)
%     example:
%           plot_itc_avg_spectrum(eeg_data,200,[0.5,4.5],1:4);
%
% Nai Ding & Wen Zhang, 2016
% [e-mail]
 
f=1:size(eeg_data,1);f=f-1;f=f/size(eeg_data,1);f=f*fs;
az = zeros(size(eeg_data,1),size(eeg_data,3),size(eeg_data,2));
for ch = 1:size(eeg_data,2)
  az(:,:,ch) = angle(fft(squeeze(eeg_data(:,ch,:))));
end
itc_sprectrum = squeeze(pcoh3(az));
avg_spectrum = abs(fft(mean(eeg_data,3)));
 figure;
 subplot(211);
 plot(f,itc_sprectrum);
 xlim(display_frequency_range);
 xlabel('frequency (Hz)')
 ylabel('inter-trial phase coherence')
 try
  set(gca,'xtick',f_label);end
 title('itc\_sprectrum');
 subplot(212);
 plot(f,avg_spectrum);
 xlim(display_frequency_range);
 title('avg\_spectrum');
 xlabel('frequency (Hz)')
 ylabel('power (a.u.)')
 try
  set(gca,'xtick',f_label);end
end
 
function r=pcoh3(ag)
c=cos(ag);
s=sin(ag);
r=mean(c,2).^2+mean(s,2).^2;
end

When I run plot_itc_avg_spectrum on its own, with a MatLab array of the form time x channel x trials as the input, it works fine, and gives me a plot showing a frequency spectrum of the data with correctly labelled axes.

However, when I run the whole seedssE script, with various functions designed to perform Denoising Source Separation (DSS), the script runs without any obvious errors, but produces a blank plot.

Both of these plots are attached.

How can I ensure that the plot that results from running the whole seedssE script is not blank, and looks like the plot resulting from running plot_itc_avg_spectrum on its own?

[Star Strider]

The most common reasons for a blank plot, in my experience at least, are:

(1.) All the data are NaN, since one NaN value in a vector input to the fft function will result in the entire vector being NaN. The NaN values may not be in the original data (and likely are not, since it plots correctly otherwise), although could be introduced during the additional processing.

(2.) Plotting points in a loop without using a linestyle with a marker, since MATLAB requires a line segment of at least two points be defined to draw a lline between them.

We have no way to determine if these or other possible problems could be causing the blank plots.

.

[Alicia Mason]

Thank you for your response!

I’ve had a close look at the script, and I think I have some idea as to where it might be going wrong.

In the function for plotting the frequency spectrum, plot_itc_avg_spectrum, the input should be MEG data in the form of a 3D array of the shape time*channel*trial.

However, when this function is called in the seedssE script, the input to the function is the variable y, which is created as follows:

y=fold(unfold(xf)*todss(:,1:10)*fromdss(1:10,:), size(xf,1));

xf is an array of size 12800x194x2, todss is a matrix of size 194x161, fromdss is a matrix of size 161x194, and unfold() is a function that will turn a 3D array into a 2D matrix - in the case of xf, it goes from 12800x194x2 to 25600x194.

I looked up the definition of fold() in the MatLab documentation, and it’s complicated, but it seems that the first argument should be a function that is called on the first two elements of the second argument (a vector) and then repeatedly called on the result and the next until all elements of the vector are completely combined.

In this call to fold(), however, the first argument unfold(xf)*todss(:,1:10)*fromdss(1:10,:) is not a function, and size(xf,1) is simply a 1 by 1 vector with the value 12800. Indeed, this is the value that ends up being assigned to the variable y, which is also a 1 by 1 vector of value 12800. Therefore, this is not the correct input to plot_itc_avg_spectrum, which requires a 3D array of shape time*channel*trial.

My big question now is the following: would this fold() function work if it were provided with a different input, or will it always produce a 1 by 1 vector regardless of the input? If the latter is true, how can I change the code so that I provide the correct input to plot_itc_avg_spectrum? I have tried using the variable xf, since this is of the correct shape, and while this does plot a readable graph, this variable has not undergone the denoising processes elsewhere in the script, and so the data displayed is still noisy.

Do you have any thoughts about how best to resolve this problem?

Thank you!

[Star Strider]

My pleasure!

The problem of course is not only the symbolic fold function (that I doubt is appropriate here since the arguments do not appear to be symbolic) but also ‘unfold’ that is not a MATLAB function, so I believe these are actually both proprietary functions that are not provided here.

I have no idea what the code does, or what the problem could be. I suggested two possibilities earlier, and those remain the only possibilities I can think of here.

You likely have to get inside the code and see if it is creating NaN values, since I doubt that it is plotting in a loop (although that remains a possibility).

Also consider writing the authors to see if they can solve this. Yours would not be the first (or only) time that published methods are not reproducible.

.

[Alicia Mason]

It seems like the fold() function was indeed a proprietary function, just like unfold() - I managed to track this down, and now it works! Thank you so much for giving me the clue to this!

[Star Strider]

My pleasure!

I still have no idea what ‘fold’ and ‘unfold’ do in that code, or even what they are.

.