Data and Code spatiotemporal white noise

spatiotemporalWhiteNoise/code spatiotemporal white noise/BC_NL.m

@@ -0,0 +1,73 @@
+function [NL_output]=BC_NL(STA_output)
+%[NL_output]=BC_NL(STA_output)
+%this function computes the output nonlinearity (NL) from the linear-nonlinear model (LN-model)
+%for continous voltage signals of bipolar cells (BC).
+%Inputs:
+%               STA_output = output of the BC_STA function; containing the
+%               STA and stimulus signal
+%
+%Outpus:        NL_output = structure containing the output nonlinearity
+%
+%
+% Note: the function is built only for cells recorded without frozen frames.
+% But additional comments are made for cells with frozen frames.
+% code by HS, last modified Feb 2021
+
+%-------------------------------------------------------------------------------
+%% 1. convolve zoom STA with stimulus (to get the generator signal)
+%get the linear filter (normalized, see Methods of Schreyer&Gollisch,2021)
+linearFiler=normSTA(STA_output.STA2D_zoom);
+%get stimulus signal
+stimulus=STA_output.stimulus2D_zoom; 
+%%if you have frozen frames, only the stimulus signal of the running frames
+%%is used for the convolution. Also the stimulus signal has to be convolved
+%%for each trial of the running noise separately (because it is not a
+%%continous signal -> it was interupted by the frozen noise).
+
+%do the convolution with the stimulus signal
+generator_signal=filter2(linearFiler,stimulus,'valid');
+%%IMPORTANT: if you have a best spatial and temporal component as described in Point 7
+%in the function BC_STA.m, this part hast to be done in two lines, e.g.:
+%NL_PredXaxis2=filter2(sp_bestnorm,stimulus,'valid'); %first
+%convolve the spatial dimension
+%generator_signal=filter2(temp_bestnorm,NL_PredXaxis2,'valid'); %then the
+%temporal dimension
+%%Also note, to avoid noise contributions from pixels outside the receptive
+%%field,it is important to cut the STA into a smaller region around the maximum pixel, this
+%%gives a much cleaner NL-> see BC_STA.m and Methods in Schreyer&Gollisch (2021).
+    
+%-------------------------------------------------------------------------------
+%% 2. get the y axis of the NL
+%get the voltage responses
+y_axisNL_unsorted=STA_output.membranPotAVG(STA_output.STA_fnbr:end); %sthe first value is at the length of the filter
+
+%-------------------------------------------------------------------------------
+%% 3. sort the axis and do the binning
+%sort the x axis
+[NL_xAxis,indx]=sort(generator_signal);
+%sort the corresponding y axis
+NL_yAxis=y_axisNL_unsorted(indx);
+%do a binning with percentile to have same amount of data in each bin
+%make 40 bins
+nbrbins=40;
+bins= doBin(NL_xAxis,NL_yAxis,nbrbins);
+
+%-------------------------------------------------------------------------------
+%% 4. plot the nonlinearity
+figure;
+plot(NL_xAxis,NL_yAxis,'.','color','k') %plot non-binned signal
+hold on;
+plot(bins.NL_xbin,bins.NL_ybin,'o','markerfacecolor','r','markeredgecolor','r') %plot binned signal
+title(['NL: filename: ' STA_output.filename])
+axis tight
+xlabel('generator signal')
+ylabel('voltage (mV)')
+
+%-------------------------------------------------------------------------------
+%% 5. add variables into the output
+NL_output.NL_xAxis=NL_xAxis;
+% STA_output.STA=STA; %not  normalized
+NL_output.NL_yAxis=NL_yAxis;
+NL_output.bins=bins; %for reshapeing to 3D structure if needed
+
+end

spatiotemporalWhiteNoise/code spatiotemporal white noise/BC_STA.m

@@ -0,0 +1,145 @@
+function [STA_output]=BC_STA(loadData,filename,pixelSize,seed,Hz)
+%[STA_output]=BC_STA(loadData,filename,pixelSize,seed)
+%this function loads the voltage trace of the selected bipolar cell (BC)
+%and perorms a response-weighted average analysis analogous to the commonn
+%spike-triggered average (STA). For simplicity we refer to the analysis
+%STA,even though we work with a continous voltage signal and not spikes. 
+%Inputs:        
+%               loadData = the folder path of the white noise files
+%               filename = the filename of the BC of interest
+%               pixelSize = the size of a square (checker) in monitor
+%               pixels (see Documentation_Data.pdf)
+%               seed = seed for the random number generator
+%               Hz = stimulus update rate
+%
+%Outpus:        STA_output = structure containing the STA, stimulus signal,
+%               averaged membrane potential for calculating the output
+%               nonlinearity
+%
+%
+% Note: the function is built only for cells recorded without frozen frames.
+% But additional comments are made for cells with frozen frames.
+% code by HS, last modified Feb 2021
+
+%% constants
+%OLED dimensions for regenerating the stimulus
+monitorXPixels=800; %x pixels of OLED monitor; 
+monitorYPixels=600; %y pixels of OLED monitor; 
+microPerPix=2.5; %one pixel=2.5 micrometer in size
+STA_TimeBack=1000; %time in ms you want to go back for generating the STA (to get smaller data I use here 1000ms, in the manuscript we use 2000ms). 
+
+%-------------------------------------------------------------------------------
+%% 1.load cell
+data=load([loadData,filename]); %load the file
+
+%-------------------------------------------------------------------------------
+%% 2. regenerate the stimulus
+Nframes=length(data.ttls_ms); %number of frames shown (how many times the stimulus was updated) 
+Nx=ceil(monitorXPixels/pixelSize); %number of squares shown X dimension
+Ny=ceil(monitorYPixels/pixelSize); %number of squares shown Y dimension
+nbrCheckers_total=Nx*Ny; %total number of checkers/square presented on the screen
+stimulus = recreateBinaryWhiteNoiseStimulus(Nx, Ny, Nframes, seed); %the function is available under: https://github.com/gollischlab/RecreateWhiteNoiseStimuliWithMatlab 
+stimulus =stimulus(:,:,1:end-1);%take the last frame out ->needed for doing the average per frames see avgPerFrame
+%the output stimulus has 3 dimensions: Y,X (nbr of checkers on Y and X axis) and time. Easiest to understand the stimulus is by
+%plotting for example the first stimulus frame presented on the retina:
+%imagesc(stimulus(:,:,1));
+%Note: the stimulus values are drawn onto the
+%retina from the lower left corner, yet matlab displays it from the upper left corner.
+%Therefore the y axis has to be flipped when displaying.
+%set(gca,'YDir','normal');
+%For me the easiest is to directly flip the stimulus y dimension (and then I dont flip the axis when displaying):
+stimulus_origD = flipdim(stimulus,1);
+clear stimulus
+%Yet that might be personal preferences.Note: For generating the STA it does not
+%matter whether it is flipped or not. However if the STA is convolved with the natural movies it matters.
+%Further, if you have frozen and running frames, regenerate the stimulus separatly
+%for running and frozen frames (see Documentation_Data.pdf for the corresponding seeds).
+
+%-------------------------------------------------------------------------------
+%% 3. make average per stimulus frame
+membranPotAVG=avgPerFrame(data.membPot,data.ttls_ms); % weights needed for computing the STA
+ 
+%-------------------------------------------------------------------------------
+%% 4. Computing the STA
+%decide for numbers of frames you want to go back in time
+STA_fnbr=(ceil(STA_TimeBack/(1000/Hz))); % define the number of frame you want to go back in time
+%here I prefer to work with a 2D signal of the stimulus
+stimulus_2D = reshape(stimulus_origD, nbrCheckers_total, []);
+tic
+[~,STA_2D]=buildSTA ...
+    (stimulus_2D,membranPotAVG,STA_fnbr);
+toc
+%if you like the STA as a 3D e.g.:
+STA_3D=reshape(STA_2D,Ny,Nx, ...
+    size(STA_2D,2));% Y (monitor axis), X (monitor axis) and Z (STA_fnbr)
+
+%-------------------------------------------------------------------------------
+%% 5. Make a zoom-in of the STA
+%cut a window of around 100 micrometer on each side from the maximum. The
+%zoom-in of the STA is better to avoid overfitting for the NL and
+%prediction.
+oneCheckInMicron=pixelSize*microPerPix;
+zoom_RF=ceil(100/oneCheckInMicron);
+
+%get the maximum of the STA to make a zoom-in
+[max_Color,indx_Zoom]=max(abs(STA_3D(:)));
+[posY,posX,posZ]=ind2sub(size(STA_3D),indx_Zoom);
+%cut a window around the STA ->make a zoom-in STA
+%zoom STA -> here you take the x-y coordinates
+rangeY_zoom = posY+(-zoom_RF:zoom_RF); rangeY_zoom = rangeY_zoom(rangeY_zoom>0 & rangeY_zoom<=size(STA_3D,1));
+rangeX_zoom = posX+(-zoom_RF:zoom_RF); rangeX_zoom = rangeX_zoom(rangeX_zoom>0 & rangeX_zoom<=size(STA_3D,2));
+STA_3D_zoom=STA_3D(rangeY_zoom,rangeX_zoom,:); %take the data from the 2000ms STA
+%reshape the STA
+STA2D_zoom=reshape(STA_3D_zoom,size(rangeY_zoom,2)*size(rangeX_zoom,2),size(STA_3D_zoom,3));
+%get the zoom in of the stimulus signal
+stimulus3D_zoom=stimulus_origD(rangeY_zoom,rangeX_zoom,:); 
+%reshap to 2D
+stimulus2D_zoom=reshape(stimulus3D_zoom,size(rangeY_zoom,2)*size(rangeX_zoom,2),size(stimulus3D_zoom,3));
+
+%-------------------------------------------------------------------------------
+%% 5. visualize STA: plot the frame with maximal intensity value
+%I also like to plot all/multiple STA frames around the maximum.
+colorLimits = [-max_Color,max_Color]; %mainly relevant if you plot multiple frames and want to keep the same colorbar
+%plot
+figure;
+imagesc(STA_3D(:,:,posZ))
+caxis(colorLimits);
+colormap('gray'); %how the sitmulus was shown
+axis equal; axis tight;
+title(['STA best frame, filename: ' filename]) 
+
+%-------------------------------------------------------------------------------
+%% 6. add variables into the output
+STA_output.STA_2D=STA_2D;
+STA_output.STA2D_zoom=STA2D_zoom;
+% STA_output.STA=STA; %not  normalized
+STA_output.stimulus2D_zoom=stimulus2D_zoom;
+STA_output.Nx=Nx; %for reshapeing to 3D structure if needed
+STA_output.Ny=Ny;
+STA_output.membranPotAVG=membranPotAVG;
+STA_output.STA_fnbr=STA_fnbr;
+STA_output.filename=filename;
+STA_output.rangeX_zoom=rangeX_zoom;
+STA_output.rangeY_zoom=rangeY_zoom;
+STA_output.pixelSize=pixelSize;
+
+%-------------------------------------------------------------------------------
+%% 7. additional comments:
+%From the zoom STA, I seperate the response into the highest-ranked spatial and temporal
+%components by singular-value decomposition (see Mansucript Schreyer & Gollisch, 2021)
+%example code for the SVD:
+%[sp,sigma,temp] = svd(STA,'econ'); %use the non-normalized STA here!
+%sp=spatial component and temp= the temporal component. 
+%the highest-ranked spatial and temporal component has to be defined e.g. by manually looking
+%at the components.
+%From here, I fitted a 2D Gaussian to the best spatial component.
+%Important for nonlinearity calculation:
+%To avoid noise contributions from pixels outside the receptive field, I reduced 
+%the number of elements in the STA by setting pixel values of the STA outside 
+%the 3-sigma contour of the Gaussian fit to zero and re-separating the spatiotemporal 
+%receptive field within this window into the highest-ranked spatial and temporal components 
+%by singular-value decomposition. Each component was then normalized to unit Euclidean norm. E.g. 
+%sp_bestnorm=normSTA(sp_best);
+%temp_bestnorm=normSTA(temp_best);
+
+end

spatiotemporalWhiteNoise/code spatiotemporal white noise/BC_mainWhiteNoise.m

@@ -0,0 +1,26 @@
+%Guide for performing a linear-nonlinear model (LN-model)
+%analysis from voltage recordings of bipolar cells (BC) under spatiotemporal white noise
+%Check out the Documentation_Data.pdf, which contains relevant information
+%for analysing the white noise stimulus
+%check also the paper from Chichilnisky, E.J. (2001). 'A simple white noise 
+%analysis of neuronal light responses', Network, 12, 199-213.
+%it explains the LN-model for spiking neurons nicely! Here we adapted this
+%method to non-spiking neurons (continous voltage signal).
+
+clear;
+close all;
+
+%% First LN-Model Component: Filter (spatiotemporal receptive field =STA)
+loadData='Z:\DataAvailable\spatiotemporalWhiteNoise\'; %path of the BC files for spatiotemporal white noise (has to be adapted to your own path)
+filename='16d200015Comp.mat'; %BC file to analyse: for example here we use file 155200029Comp from cell 2, recorded in retina 1 on 150520
+pixelSize=9; %is given in the excel file LogInfo_SpatioTemporalWhiteNoise.xlsx for the corresponding cell
+seed=-10000; %see also documention_data.pdf
+Hz=25; % given in the excel file LogInfo_SpatioTemporalWhiteNoise.xlsx for the corresponding cell
+[STA_output]=BC_STA(loadData,filename,pixelSize,seed,Hz);%example code receptive field analysis=linear filter
+
+%% Second LN-Model Component: Output nonlinearity
+[NL_output]=BC_NL(STA_output);%example code output nonlinearity analysis
+
+%% Prediction of the response with the LN-Model (filter and output nonlinearity):
+[pred_output]=BC_pred(STA_output);%example code prediction
+

spatiotemporalWhiteNoise/code spatiotemporal white noise/BC_pred.m

@@ -0,0 +1,142 @@
+function [pred_output]=BC_pred(STA_output)
+%[pred_output]=BC_pred(STA_output)
+%this function shows the prinicpal calculation for the the prediction to 
+%the white noise stimulus with the linear-nonlinear model (LN-model)for 
+%continous voltage signals of bipolar cells (BC).
+%Inputs:
+%               STA_output =  output of the BC_STA function; containing the
+%               STA and stimulus signal      
+%
+%Outpus:        pred_output = structure containing the output of the
+%               prediction
+%
+%
+% Note: the function is built only for cells recorded without frozen frames.
+% But additional comments are made for cells with frozen frames.
+% code by HS, last modified Feb 2021
+
+%-------------------------------------------------------------------------------
+%% 1. recompute STA and nonlinearity without the stimulus signal that we want to predict
+%%Important we want to predict a stimulus sequence that was not used to
+%%build the model (filter and nonlinearity). Therefore we have to redo the
+%%filter(STA) and nonlinearity without the stimulus-sequence that we want 
+%to predict the response(see Methods of Schreyer&Gollisch,2021).
+%If we recorded Frozen frames, we predict the response to the frozen
+%frames from the model (filter and nonlinearity) we build from the running
+%frames.
+
+%% 1.1 STA
+%E.g. to make it easy: we want to predict the response to the last 300frames of the white
+%noise stimulus
+%take out the stimulus sequence you want to predict
+stimPred=STA_output.stimulus2D_zoom(:,end-299:end); %the last 300frames
+%voltage trace to predict
+voltPred=STA_output.membranPotAVG(end-299:end); %the response to the last 300 frames
+%build the STA without the last 300 frames
+[STA_norm]=buildSTA ...
+    (STA_output.stimulus2D_zoom(:,1:end-300),...
+    STA_output.membranPotAVG(1:end-300),STA_output.STA_fnbr);
+STA3D_norm=reshape(STA_norm,size(STA_output.rangeY_zoom,2),size(STA_output.rangeX_zoom,2), ...
+    size(STA_norm,2));% Y-zoom checkers (monitor axis), X-zoom checkers (monitor axis) and Z (STA_fnbr)
+
+%% 1.2 Nonlinearity
+%build the NL from the data without the last 300 frames
+generator_signal=filter2(STA_norm,STA_output.stimulus2D_zoom(:,1:end-300),'valid');
+%get the corresponding membrane potential responses
+y_axisNL_unsorted=STA_output.membranPotAVG(STA_output.STA_fnbr:end-300); %the first value is at the length of the filter
+
+%sort the x axis
+[NL_xAxis,indx]=sort(generator_signal);
+%sort the corresponding y axis
+NL_yAxis=y_axisNL_unsorted(indx);
+%do a binning with percentile to have same amount of data in each bin
+%make 40 bins
+nbrbins=40;
+bins= doBin(NL_xAxis,NL_yAxis,nbrbins);
+
+%-------------------------------------------------------------------------------
+%% 2. compute the prediction
+%% 2.1 convolution
+%convolve the new stimulus sequence (not used for prediciton) with the
+%normalized STA
+stim2pred=filter2(STA_norm,stimPred,'valid');
+
+%% 2.2 compute the prediction through linear interpolation
+NL_x=bins.NL_xbin;
+NL_y=bins.NL_ybin;
+%example of how to make the prediction through linear interpolation
+predMemP_y=nan(1,size(stim2pred,2));
+for kk=1:size(stim2pred,2)
+    nbrTemp=stim2pred(kk); %first x-axis value to predict
+    [xIndxS,xIndxL]= getIndex(NL_x,nbrTemp,1); %get those values that are on left and rigth
+    xValueS=NL_x(xIndxS);
+    xValueL=NL_x(xIndxL);
+    %do the linear interpolation
+    [polyCoeff] = polyfit([xValueS;xValueL],[NL_y(xIndxS);...
+        NL_y(xIndxL)],1);
+    predMemP_y(kk)=polyCoeff(1)*nbrTemp+polyCoeff(2);
+end
+
+%-------------------------------------------------------------------------------
+%% 3. evaluate the prediction performance
+%pearson correlation
+corrcoeff_temp=corrcoef(voltPred(STA_output.STA_fnbr:end),predMemP_y); %the first value is at the length of the filter
+pearCorrSqrt=corrcoeff_temp(2)^2;
+
+%-------------------------------------------------------------------------------
+%% 4. plot the output
+figure;
+%STA
+subplot(2,2,1)
+[~,indx_Zoom]=max(abs(STA3D_norm(:))); %max intensity value of STA
+[~,~,posZ]=ind2sub(size(STA3D_norm),indx_Zoom); %get the frame position
+imagesc(STA3D_norm(:,:,posZ))
+colormap('gray'); %how the sitmulus was shown
+axis equal; axis tight;
+title('STA best frame')
+%NL
+subplot(2,2,2)
+plot(NL_xAxis,NL_yAxis,'.','color','k') %plot non-binned signal
+hold on;
+plot(bins.NL_xbin,bins.NL_ybin,'o','markerfacecolor','r','markeredgecolor','r') %plot binned signal
+title('NL')
+axis tight
+xlabel('generator signal')
+ylabel('voltage (mV)')
+%prediction
+subplot(2,2,3:4)
+plot(voltPred(STA_output.STA_fnbr:end),'k') %original membrane potential
+hold on;
+plot(predMemP_y,'r')%predicted membrane potential
+title(['pCorr:' mat2str(round(pearCorrSqrt,2)) ' Filename: ' STA_output.filename])
+xlabel('stimulus frames')
+ylabel('voltage (ms)')
+legend('original response','predicted response')
+
+%-------------------------------------------------------------------------------
+%% 5. add variables into the output
+pred_output.pearCorrSqrt=pearCorrSqrt;
+
+end
+
+function [IndS,IndL]= getIndex(x,valConv,Points)
+IndS=find(x<valConv,Points,'last'); %find the last 5 closest smaller values
+IndL=find(x>=valConv,Points);
+if isempty(IndS)|| length(IndS)<Points %if you are on the left border
+    if isempty(IndS) %take value from 1:10
+        IndL=find(x>=valConv,Points*2);
+        IndS=IndL(1);
+        IndL=IndL(2:end);
+    else
+        IndL=find(x>=valConv,Points*2-length(IndS));
+    end
+elseif isempty(IndL)|| length(IndL)<Points %if you are on the right border
+    if isempty(IndL) %take value from 1:10
+        IndS=find(x<valConv,Points*2,'last');
+        IndL=IndS(end);
+        IndS=IndS(1:end-1);
+    else
+        IndS=find(x<valConv,Points*2-length(IndL));
+    end
+end
+end

spatiotemporalWhiteNoise/code spatiotemporal white noise/avgPerFrame.m

@@ -0,0 +1,18 @@
+function [membranPotAVG]=avgPerFrame ...
+        (membranePot,ttls_ms)
+%function [membranPotAVG]=sharpEL_buildAveragePerFrame ...
+%         (membranePot,ttls_ms_correctedPulse)
+%builds the average membrane potential between frames for the white noise
+% -> same principle for frozen frames
+% INPUT:        membranePos -> membrane Potential
+%               ttls_ms -> pulses in ms
+% OUTPUT:       membranPotAVG -> average membrane potential per frames
+
+membranPotAVG=zeros(1,length(ttls_ms)-1);%preallocate
+
+%make the average membrane potential between two stimulus frames
+for jj=1:length(ttls_ms)-1 
+    membranPotAVG(jj)=mean(membranePot(ttls_ms(jj) ...
+        :(ttls_ms(jj+1)-1)));%detail: -1ms because the next ms is already the start of the new stimulus frame
+end
+end

spatiotemporalWhiteNoise/code spatiotemporal white noise/buildSTA.m

@@ -0,0 +1,34 @@
+function [STA_mean_normalized,STA_mean]=buildSTA ...
+    (stim_signal,membraneP,STA_fnbr)
+% [STA_mean_normalized,STA_mean]=buildSTA ...
+%    (stim_signal,membraneP,STA_fnbr)
+%is getting the contrast value that happen before each response
+%and multiplies them with the average membrane Potential per frame.
+% INPUT:        stim_signal -> stimulus signal containing the contrast values in 2D
+%               membraneP -> membrane potential average per frame
+%               STA_fnbr -> nbr of frames you go back in time
+% OUTPUT:       STA_mean_normalized -> normalized STA, relevant for
+%               computing the NL
+%               STA_mean -> not normalized STA
+
+%------------------------------------------------------------------------------------------
+%% 1. Loop over all membrane potential values
+STA_fnbr=STA_fnbr-1; %the first value is one frame, so to take exactly 10 frames, substract one here.
+for kk=1:length(membraneP)-(STA_fnbr)
+    if kk==1
+        signal1=membraneP(kk+(STA_fnbr))*stim_signal(:,kk:kk+(STA_fnbr)); %multiply the average values with the stimulus signal before; it start now from 1 and goes until the element kk
+    else
+        signal2=membraneP(kk+(STA_fnbr))*stim_signal(:,kk:kk+(STA_fnbr)); %multiply the average values with the stimulus signal before
+        signal1=signal1+signal2; %make the sum to get the average at the end
+        clear signal2;
+    end
+end
+%build the average
+STA_mean=signal1/length(1:length(membraneP)-(STA_fnbr)); %not normalized
+
+%normalize
+STA_mean_normalized = STA_mean(:) / sqrt(sum(STA_mean(:).*STA_mean(:))); %relevant for computing the NL
+STA_mean_normalized=reshape(STA_mean_normalized,size(STA_mean));
+
+
+end

spatiotemporalWhiteNoise/code spatiotemporal white noise/doBin.m

@@ -0,0 +1,34 @@
+function  [bin]=doBin(NL_x,NL_y,xbin)
+%this subfunction is doing the binning for nonLinearity
+%get different bins
+
+%% 1. get equal amount of observation inside a bin
+%standard score:
+x=prctile(NL_x,linspace(0,100,xbin+1)); %you want equal amount of observations in each bin!And not one bin of only e.g. 10 observations.
+%I want to now where the percentil are
+%prctile gives you now where the edges of the corresponding percentile rank
+%are.
+
+x(end)=x(end)+x(end)/100; %the last bin should be a bit bigger so that the last observation
+%can also go inside it
+[observations,indx]=histc(NL_x,x);
+
+%little trick so that you have 40 equal bin and not one bin at the end
+%that is empty or just has one value
+observations=observations(1:end-1);
+%make average in y axis
+NL_ytemp=cell(1,length(observations));
+for mm=1:length(NL_ytemp)
+    NL_ytemp{mm}=NL_y(indx==mm);
+end
+%make average in x axis
+x_AVG=cell(1,length(observations));
+for mm=1:length(x_AVG)
+    x_AVG{mm}=NL_x(indx==mm);
+end
+
+%get the average values per bin
+bin.NL_ybin=cellfun(@nanmean,NL_ytemp);
+bin.NL_xbin=cellfun(@nanmean,x_AVG);
+
+end

spatiotemporalWhiteNoise/code spatiotemporal white noise/normSTA.m

@@ -0,0 +1,5 @@
+function normX =normSTA(x)
+%this function normalizes (with euqlidien distance) and reshapes a vector to its original size
+normX = x(:) / sqrt(sum(x(:).*x(:)));
+normX=reshape(normX,size(x));
+end