The following is the code:-
clear all;
close all;
clearvars;
tic
%% Read the data
load data2mw
figure();
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
time=data(:,1);
voltage=data(:,2);
t1=2.1755
t2=2.1770
i1=find(time(:,1)>t1,1)
i2=find(time(:,1)>t2,1)
Time=time(i1:i2);
Volburst=voltage(i1:i2);
figure()
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');
figure()
plot(Time,Denoised_Burst_signal)
%% calculate the analytical signal and get the envelope %%
analytic = hilbert(Denoised_Burst_signal);
env = abs(analytic);
figure(10);plot(env)
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);
figure(11)
plot(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));
h=1;
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));
h=h+1;
else
% ----------- update noise --------------%
if mean(env(i:i+k)) < THR_SIG
nois = mean(env(i:i+k));
else
if ~isempty(nois_buf)
nois = mean(nois_buf);
end
end
end
thres_buf(i) = threshold;
nois_buf(i) = nois;
% ------------------------- update threshold -------------------------- %
if h > 1
THR_SIG = nois + 0.50*(abs(threshold - nois));
end
THR_buf(i) = THR_SIG;
end
figure,ax(1)=subplot(211);
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',...
'orientation','horizontal');
grid on;axis tight;
ax(2)=subplot(212);plot(env);
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');
linkaxes(ax,'x');
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);
%%
x=1:1:N;
for i = 1:length(x)
%f_est_1 = zeros(1,N);
yy = Volburst(digitalCellNew{1,i})
tt = Time(digitalCellNew{1,i})
figure()
plot(tt,yy)
%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
%N=nfft;
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
figure()
plot(f_scale(1:length(y3))',y3);
f_est_1(i)= f_scale(k); % dominant frequency estimate for each burst
hold on; % Prevent image from being blown away
plot(f_est_1,v,'ro');
axis('auto'),grid('on');
xlim([0 2e+07]);
title(''),xlabel('Frequency (Hz)'),ylabel('Voltage (V)');
%
end
%%
ElaspedTime1 = toc