Converting for loops to parfor format gives error "Error: Unable to classify the variable 'f_est_1' in the body of the parfor-loop."

2 views (last 30 days)
Hello, I am trying to convert the for-loop into parfor loop for faster computation (because of huge amount of data).
The code and data (sample data) is attached for the reference. The code works fine when compiled with for-loop (attached code is with for-loop).
When I change the for-loop to parfor I get the error "Error: Unable to classify the variable 'f_est_1' in the body of the parfor-loop".
The attached code segments the signal based on the adaptive thresholding. Then plots the segmented data and performs the fft of it and stores the dominant frequency from each fft plots in a vector. These tasks are done using the "for-loop". I wanted to convert the for-loops to parfor, I tried changing the the way variables to be used as per the parfor variable classification but I am not bale to crack it.
how to perform this task using parfor? any help is highly appreciable.
  1 Comment
Rajendra  Rajak
Rajendra Rajak on 15 Jul 2021
The following is the code:-
clear all;
close all;
%% Read the data
load data2mw
plot(data(:,1), data(:,2)-mean(data(:,2)));
xlabel('Time (s)','fontsize', 10); ylabel('Signal magnitude (V)', 'fontsize', 10);
title('Raw Data signal', 'fontsize', 10);set(gca,'FontSize',10);
%% Extract the data for analysis
plot (Time,Volburst-mean(Volburst));xlabel('Time (seconds)');ylabel('Voltage (V)');
%% Denoise the signal
Denoised_Burst_signal = wdenoise(Volburst,12, ...
'Wavelet', 'sym8', ...
'DenoisingMethod', 'UniversalThreshold', ...
'ThresholdRule', 'Soft', ...
'NoiseEstimate', 'LevelIndependent');
%% calculate the analytical signal and get the envelope %%
analytic = hilbert(Denoised_Burst_signal);
env = abs(analytic);
k = 20;
%% moving average of the analytical signal
Smooth_window = 40;
env = conv(env,ones(1,Smooth_window)/Smooth_window);
env = env(:) - mean(env);
env = env/max(env);
%% Threshold and initialize the buffers
THR_SIG = 4*mean(env);
nois = mean(env)*(1/3);
threshold = mean(env);
thres_buf = zeros(1,length(env)-k);
nois_buf = zeros(1,length(env)-k);
THR_buf = zeros(1,length(env));
LogicalSignal =zeros(1,length(env));
for i = 1:length(env)-k
%------ update threshold 10% of the maximum peaks found ------------%
if env(i:i+k) > THR_SIG
LogicalSignal(i) = max(env);
threshold = 0.1*mean(env(i:i+k));
% ----------- update noise --------------%
if mean(env(i:i+k)) < THR_SIG
nois = mean(env(i:i+k));
if ~isempty(nois_buf)
nois = mean(nois_buf);
thres_buf(i) = threshold;
nois_buf(i) = nois;
% ------------------------- update threshold -------------------------- %
if h > 1
THR_SIG = nois + 0.50*(abs(threshold - nois));
THR_buf(i) = THR_SIG;
plot(Volburst),hold on,plot(LogicalSignal.*max(Volburst),'r','LineWidth',0.5),
hold on,plot(THR_buf,'--g','LineWidth',0.5);
title('Raw Signal and detected Onsets of activity');
legend('Raw Signal','Detected Activity in Signal',...
'Adaptive Treshold',...
grid on;axis tight;
hold on,plot(THR_buf,'--g','LineWidth',0.5),
hold on,plot(thres_buf,'--r','LineWidth',0.5),
hold on,plot(nois_buf,'--k','LineWidth',0.5),
title('Smoothed Envelope of the signal(Hilbert Transform)');
legend('Smoothed Envelope of the signal(Hilbert Transform)',...
'Adaptive Treshold',...
'Activity level','Noise Level','orientation','horizontal');
zoom on;
axis tight;
grid on;
%% Calcuation of burst parameters and velocity prediction
digitalCell = bwconncomp(LogicalSignal);
%digitalCell.num = cellfun(@length,digitalCell.PixelIdxList)
CellLength = cellfun(@length,digitalCell.PixelIdxList);
CellLength2keep = CellLength>150;
N= sum(CellLength2keep);
digitalCellNew = digitalCell.PixelIdxList(CellLength2keep);
for i = 1:length(x)
%f_est_1 = zeros(1,N);
yy = Volburst(digitalCellNew{1,i})
tt = Time(digitalCellNew{1,i})
%format long
% FFT of the each burst
% FFT of the Burst Voltage
L=length(tt); % Length of signal
nfft=2^nextpow2(L); %calculates the nearest power of 2 w.r.t. l
Ts=tt(end)-tt(1); % Total Observation time
Fs=L/Ts; % Sampling frequency (will be same as fps)
% F_nyquist=Fs/2; % Nyquist frequency
yy_in_PRIME = yy - mean(yy);
y1 = fft(yy_in_PRIME,L); % Discrete Fourier transform (DFT)
y2 = 2*((abs(y1))/L); % Amplitude of the DFT
[v,k] = max(y2); % find maximum value (v) and it's position in the array (k)
% f_totalscale = (0:(N-1))* Fs/N; % total frequency scale
f_scale=(0:(L-1)/2)* (Fs/L); % Half the total frequency range
y3=y2(1:round((L/2))); % half the power spectrum
f_est_1(i)= f_scale(k); % dominant frequency estimate for each burst
hold on; % Prevent image from being blown away
xlim([0 2e+07]);
title(''),xlabel('Frequency (Hz)'),ylabel('Voltage (V)');
ElaspedTime1 = toc

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 15 Jul 2021
f_est_1(i)= f_scale(k); % dominant frequency estimate for each burst
Inside parfor i, at any one time, a (possibly discontinuous) subset of f_est_1 has had valid values written to it. parfor does not execute in sequential order (deliberately), so parfor will not write to f_est_1(1) then f_est_1(2) and so on. It might happen to write to f_est_1(17) then f_est_1(11) then f_est_1(20) then f_est_1(3) and so on. All you know about the order is that by the end of the loop, all of the entries have been filled in.
hold on; % Prevent image from being blown away
But there you try to plot using all of f_est_1 including the parts that have not been filled in yet.
Rajendra  Rajak
Rajendra Rajak on 19 Jul 2021
Hi Walter,
parfor i = 1:length(LL)
%------ update threshold 10% of the maximum peaks found ------------%
if env(i:i+k) > THR_SIG
LogicalSignal(i) = max(env);
thres_buf(i) = threshold;
% ------------------------- update threshold -------------------------- %
THR_buf(i) = THR_SIG;
This part of the code consumes more time and also it gives warning " the entire array or struct 'env' is a broadcast variable. This might result in unnecessary communication overhead".
Is there a way to get rid of this? (for 12 Mega data points it took 85 minutes only for this part)

Sign in to comment.

More Answers (1)

Raymond Norris
Raymond Norris on 15 Jul 2021
I'm gathering the for-loop you're trying to re-write is
for i = 1:length(x)
Assuming you're able to re-write this, have you considered what will happen to plotting inside the parfor? Keep in mind, the entire body of code in the parfor is running elsewhere (i.e. the workers). Therefore, you won't see plotting. You have two options for plotting "parallelization" (I put this in quotes since it will need to happen outside of the parfor, where you are running serial code).
  • DataQueue: look at creating a data queue to send data from the workers to the client, which then calls a function to update your plot(s).
  • parfeval: you'll still need a parallel pool, but rather then using the built-in parfor, create tasks with parfeval and then update the plot(s) when the tasks come back.
This doesn't answer your parfor question, but if you implement parfeval, (A) you ought to be able to plot your data and (B) it'll probably solve your original issue.
Rajendra  Rajak
Rajendra Rajak on 16 Jul 2021
Thank you Walter for the descriptive reposne to solution of the problem and sharing the link.
I am going to restructure the code and get back to you in the same thread to list out any discrepancies MATALB throws during compilation.

Sign in to comment.




Community Treasure Hunt

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

Start Hunting!