INTENSITY WAVEFORM RECONSTRUCTION FROM INTENSITY

AUTOCORRELATION

FIELD OF THE INVENTION

The present invention relates generally to methods and systems for signal measurement and estimation, and specifically to techniques for measuring characteristics of very short pulses.

BACKGROUND

Intensity autocorrelation is commonly used to estimate the duration of ultra-short laser pulses. Such pulses cannot be easily measured directly by optoelectronic methods, since the response time of common detectors and measuring instruments are at best on the order of 1 ps, whereas the laser pulses themselves can be made as short as a few femtoseconds. The term "ultra-short" is used herein to refer to pulses whose duration is less than half the response time of available detectors, for example, 100 fs or less in the case of optical pulses.

The autocorrelation trace of a pulse of intensity profile I(t) is given by /^c ⁼ J ^" /(t) /(t— r)dt. The trace is usually generated by measuring nonlinear interaction between the laser pulse itself and a replica that is delayed by the variable time τ. The pulse width of I(t) can then be estimated from the width of IAC ( ^T ) > but techniques known in the art are unable to recover the actual waveform of I(t) in this manner. In fact, an infinite number of different waveforms can give rise to the same autocorrelation trace.

The present patent application relates to reconstruction of the waveform I(t). "Reconstructing a waveform," in the context of the present description and in the claims, refers to extracting the variation of the intensity I(t) as a function of time with a temporal resolution that is substantially finer than the width of the pulse itself. Thus, for example, the reconstructed waveform may include features of the temporal shape of the pulse whose width is half the total pulse width (in which case the temporal resolution is twice the pulse width), and even finer features in many cases.

A technique known as frequency-resolved optical gating (FROG) was developed to provide more precise measurements of ultra-short pulses. The technique was described by Trebino et al., in U.S. Patent 5,530,544, and in "Measuring ultrashort laser pulses in the time- frequency domain using frequency-resolved optical gating," Review of Scientific Instruments 68, pages 3277-3295 (1997). FROG is based on spectrally-resolved autocorrelation, in which the pulse gates itself in a nonlinear optical medium, and the gated part of the pulse is spectrally resolved as a function of delay. A two-dimensional phase-retrieval algorithm is then applied to recover the intensity and phase of the pulse waveform over time.

SUMMARY

Embodiments of the present invention that are described hereinbelow provide methods and apparatus for measurement of very short pulses.

There is therefore provided, in accordance with an embodiment of the invention, a method for measurement, which includes acquiring an autocorrelation trace of an ultra-short pulse of radiation. A basis of functions is defined in which an intensity profile of the pulse can be sparsely represented. The basis is searched in order to compute a sparse representation of the pulse in terms of the functions having an autocorrelation matching the autocorrelation trace. A waveform of the pulse is reconstructed using the computed sparse representation.

In some embodiments, acquiring the autocorrelation trace includes measuring an autocorrelation of a laser pulse having a duration no greater than 100 fs.

In the disclosed embodiments, the defined basis includes a set of Gauss-Hermite functions. Each of the Gauss-Hermite functions is characterized by a polynomial order, a pulse- width parameter, and a pulse-center parameter, and defining the basis typically includes selecting respective ranges of the polynomial order and the pulse-width and pulse-center parameters, and choosing the Gauss-Hermite functions within the respective ranges.

Additionally or alternatively, searching over the basis includes finding an optimal set of the functions in the basis to represent the waveform subject to a sparseness constraint. In some embodiments, the optimal set has a sparsity no greater than fifty and represents the waveform with a square error no greater than 2%. In a disclosed embodiment, the optimal set has a sparsity no greater than seventeen.

In some embodiments, reconstructing the waveform includes reconstructing features of a temporal shape of the pulse having a feature width that is less than half of a total width of the pulse.

There is also provided, in accordance with an embodiment of the invention, measurement apparatus, including an optical correlator, which is configured to measure an autocorrelation trace of an ultra-short pulse of radiation. A waveform analyzer is configured to accept a definition of a basis of functions in which an intensity profile of the pulse can be sparsely represented, to search over the basis in order to compute a sparse representation of the pulse in terms of the functions having an autocorrelation matching the autocorrelation trace, and to reconstruct a waveform of the pulse using the computed sparse representation. There is additionally provided, in accordance with an embodiment of the invention, a computer software product, including a computer-readable medium in which program instructions are stored, which instructions, when read by a processor, cause the processor to acquire an autocorrelation trace of an ultra-short pulse of radiation, to accept a definition of a basis of functions in which an intensity profile of the pulse can be sparsely represented, to search over the basis in order to compute a sparse representation of the pulse in terms of the functions having an autocorrelation matching the autocorrelation trace, and to reconstruct a waveform of the pulse using the computed sparse representation.

The present invention will be more fully understood from the following detailed description of the embodiments thereof, taken together with the drawings in which:

BRIEF DESCRIPTION OF THE DRAWINGS

Fig. 1 is a block diagram that schematically illustrates a system for pulse waveform reconstruction, in accordance with an embodiment of the present invention;

Fig. 2 is a plot that schematically represents the frequency spectrum of an ultra-short pulse;

Fig. 3 is a plot that schematically illustrates the autocorrelation trace of a pulse whose spectrum is shown in Fig. 2;

Fig. 4 is a flow chart that schematically illustrates a method for pulse waveform reconstruction, in accordance with an embodiment of the present invention;

Fig. 5 is a plot that schematically illustrates the waveform of a pulse following reconstruction by the method of Fig. 4, in accordance with an embodiment of the present invention; and

Fig. 6 is a plot showing cumulative distribution functions of reconstructed pulses as a function of sparsity level, for pulses of several different autocorrelation support durations, in accordance with an embodiment of the present invention.

DETAILED DESCRIPTION OF EMBODIMENTS

On its face, the problem of extracting the time-dependent waveform of an ultra-short pulse (or any pulse) from its one-dimensional autocorrelation trace appears to be intractable, having an infinite number of possible solutions. For this reason, approaches such as FROG collect additional information in order to transform the problem into two dimensions. Such approaches, however, entail more costly setups for measurement, followed by computation- intensive processing. Embodiments of the present invention that are described herein overcome these limitations and enable the detection of intensity waveform of an ultra-short pulse to be reconstructed and characterized solely from its one-dimensional autocorrelation trace. The waveform is reconstructed by finding a suitable functional representation of the pulse in either the time or the frequency domain. In the specific case of the autocorrelation trace, reconstruction of the waveform involves retrieval of the spectral phase of the Fourier transform

/(ω) of the intensity I(t), since the autocorrelation itself corresponds to |/(ω) | .

The present embodiments are based on the realization that ultra-short laser pulses are not arbitrary functions, but rather have a particular spectral structure. Specifically, the spectral phase of the electric field of the pulses can be described by a low-order polynomial. This characteristic of the pulse signal, taken together with the fact that the intensity of the pulses is, by definition, non-negative, is applied in embodiments of the present invention in resolving the inherent ambiguity of the autocorrelation trace.

Given the low-order polynomial spectral phase of the pulses, embodiments of the present invention are able to reconstruct the pulse waveform accurately in the time domain using a sparse set of functions within an appropriate basis, such as a set of Gauss-Hermite functions based on the Hermite polynomials. The set of functions used in this context is considered to be sparse in the sense that the waveform is represented, with an accuracy limited only by the signal/noise ratio of the autocorrelation measurement, by a sum of no more than fifty functions from the basis with non-zero coefficients. In nearly all practical cases, no more than seventeen functions from the basis are required for this purpose, and in most cases, less than ten functions are sufficient.

The embodiments of the present invention that are described hereinbelow provide a method and apparatus for measurement of ultra-short pulses that starts with acquisition of an autocorrelation trace of an ultra-short pulse in a train of ultra-short pulses. A basis of functions, such as a set of Gauss-Hermite functions, is defined in which the intensity waveform of the pulse can be sparsely represented. A search is performed over the basis in order to compute a sparse representation of the pulse in terms of the functions in the basis, such that the representation has an autocorrelation that matches the acquired autocorrelation trace. The search seeks an optimal coefficients for a set of the basis functions, subject to a sparseness constraint. The pulse waveform can be represented and output on the basis of the sparse representation that is thus computed. Fig. 1 is a block diagram that schematically illustrates a system 20 for pulse waveform reconstruction, in accordance with an embodiment of the present invention. System 20 is used to characterize pulses emitted by a mode-locked laser 22, such as a Ti: sapphire laser, which emits a train of ultra-short pulses having a consistent waveform from pulse to pulse and a pulse duration (FWHM) less than 100 fs. In the embodiment described below, the pulse duration is in the range of 30-50 fs.

For this purpose, an optical correlator 24 generates an autocorrelation trace of the pulses as follows: A beamsplitter 26 splits each pulse into two beams 28 and 30. A variable delay line 32 introduces an adjustable delay in beam 30, for example by operating an adjustable mount 36 to vary the distance between a pair of mirrors 34 from which beam 30 is reflected. Beam 28 is turned by a mirror 38 toward a lens 40, which focuses both of beams 28 and 30 into a nonlinear medium 42, such as a suitable crystal. Interaction of beams 28 and 30 in medium 42 generates higher-order radiation (typically second-harmonic radiation), which is focused by a collection lens 44 onto a high-speed detector 46, such as a suitable avalanche photodiode or photomultiplier. Alternatively, any other suitable method that is known in the art may be applied in order to measure the autocorrelation trace.

In most cases, an autocorrelation trace is measured using different pulses for different time delays. Alternatively, it is possible to measure a whole autocorrelation trace with a single pulse of sufficiently high pulse energy. In a single-shot autocorrelator, the laser beam is focused into the nonlinear crystal by a cylindrical (rather than spherical) lens, and the higher- order signal is recorded by a camera. Different spatial positions in the crystal then correspond to different time delays. For example, the Single-Shot Autocorrelator (SSA) produced by Coherent Inc. (Santa Clara, California) operates on this sort of principle.

A waveform analyzer 48 acquires and analyzes the autocorrelation trace provided by correlator 24. Analyzer 48 comprises a front-end pulse processor 50, as is known in the art, which measures the amplitudes of the electrical pulses that are output by detector 46 at each delay value (τ) set by delay line 32. A digital processor 52 records the output pulse amplitudes as a function of the delay, thus generating the autocorrelation trace /^c °f the laser pulses, and stores the trace in a memory 54 for analysis. Typically, processor 52 comprises a general purpose central processing unit (CPU), which performs this and other functions that are described herein under the control of software. The software may be downloaded to analyzer 48 in electronic form, over a network, for example, as well as stored on tangible, non-transitory computer-readable media, such as optical, magnetic, or electronic memory media. Alternatively or additionally, at least some of the functions of processor 52 may be performed by hard-wired or programmable logic circuits.

Processor 52 computes a sparse representation of the waveform of the pulses from laser 22 using a selected basis of functions, such as the set of Gauss-Hermite functions, as described in detail hereinbelow. The processor thus finds a representation of a pulse waveform having an autocorrelation that matches the autocorrelation trace, while satisfying a sparseness criterion, and reconstructs the waveform using the sparse representation that it has computed. Processor 52 typically outputs the pulse waveform and/or other characteristics via an output interface 56, such as a display or communication interface.

Fig. 2 is a plot that schematically represents the frequency spectrum, in Peta-Hertz

(PHz), of a typical ultra-short pulse of the sort emitted by laser 22. The amplitude scale is arbitrary. The pulses from the laser have a transform-limited duration of approximately 30 fs, with an added spectral chirp in this example due to passage through a glass plate 10 mm thick. The inventors have found that the spectral chirp of such pulses can be approximated well by a fifth-order polynomial expression.

Fig. 3 is a plot that schematically illustrates the autocorrelation trace of a pulse whose spectrum is shown in Fig. 2, as acquired by system 20. The trace, IAC ( ^T ) > presents the magnitude of the second-harmonic signal measured by detector 46 as a function of time delay τ between beams 28 and 30, in femtoseconds. The support of the trace in this case (i.e., the time range over which /^c is non-zero) is about 600 fs. "Non-zero" in this context means that the measured autocorrelation is significantly greater than the system noise limit, which in this case is about 0.1% of the peak autocorrelation value.

Fig. 4 is a flow chart that schematically illustrates a method for pulse waveform reconstruction, in accordance with an embodiment of the present invention. The method begins with measurement of an autocorrelation trace, such as the autocorrelation of the pulses from laser 22, as described above, at a measurement step 60. The present method is described, for the sake of clarity and concreteness, with reference to the elements of system 20 (Fig. 1) and the pulse characteristics shown in Figs. 2 and 3. The method is by no means limited to this particular context, however, and may similarly be implemented in other sorts of autocorrelation-based systems and applied to other sorts of ultra-short pulses.

In order to reconstruct the waveform of the pulses in system 20, an appropriate basis set is chosen to model the waveform, at a basis selection step 62. Given the polynomial character of the spectral chirp that is shown in Fig. 2, the inventors have found that the waveform can be represented sparsely using a basis that is built on a set of low-order polynomials in Gauss- Hermite functions. Specifically, an over-complete set of Gauss-Hermite functions can be used for this purpose. These functions are based on the Hermite polynomials, as defined by the following expression for the n-th order polynomial:

^■ d ^!l -> d Y ^l

H x) = MyVT-^-e^ = 2;r - - ) · 1

The Gauss-Hermite functions Wm _,q ^use d to model the pulses are chosen from the set: lpn _, m _,q = H _n (t)exp(-(t - t _m ) ² /M _q )

Here At _q is a pulse- width parameter (q = 1, 2, Q), t _m is a pulse-center parameter (m = 1, 2, ..., M), and the order of the Hermite polynomials is designated by n = 0, 1, 2, ..., N. This basis frame takes as its point of departure that the -spectral phase of the measured pulses can be represented in terms of low-order polynomials, but it otherwise does not assume any knowledge of the exact basis in which the temporal waveform is sparse.

For example, a frame of 180 functions may be selected by setting N = 17, M = 10, and Q = 1. The pulse-center parameters t _m are uniformly distributed over this interval. The pulse- width parameter is typically set as a function of the support of the measured autocorrelation traces, at a width selection step 64. In the present example, Ati is set to 288 fs, which is 0.36 (1/e) times the maximal support of the autocorrelation traces (800 fs). By way of comparison, for a transform-limited laser pulse (i.e., with flat phase profile), the pulse duration (FWHM) is 42 fs, while the autocorrelation support is 160 fs.

Alternatively, the selected pulse-width and pulse-center parameters may depend on the measured data. For example, the width parameter Ati may be set to a certain fraction of the width of the actual autocorrelation trace (Fig. 3), for example, 0.7 times the full width at half- maximum (FWHM) of the autocorrelation trace. Reconstruction results may be further improved (possibly at the cost of increased computational resources) by increasing the number of functions in the frame (i.e., increasing N, M and/or Q).

Analyzer 48 is programmed with the selected Gauss-Hermite basis set, and applies this set in reconstructing the pulse intensity profile, at a waveform reconstruction step 66. For this purpose, processor 52 applies a sparsity-based phase retrieval algorithm in order to find a set of coefficients of the Gauss-Hermite functions in the basis set that will reproduce the autocorrelation trace to within a target level of error. The error bound is typically set so that the square error (SE) of the reconstructed autocorrelation relative to the measured trace is no greater than a predefined limit (for example, 2%).

Specifically, at step 66, processor searches over the selected basis set to find the optimal sparse set of coefficients that will reproduce the autocorrelation trace. One algorithm that may be used for this purpose is GESPAR (GrEedy Sparse PhAse Retrieval), as described by Shechtman et al., in "GESPAR: Efficient Phase Retrieval of Sparse Signals," IEEE Transactions on Signal Processing 62, pages 928-938 (2014), which is incorporated herein by reference. GESPAR is based on a local search method with non-linear optimization of a sparsity-constrained formulation. The local-search attempts to minimize an objective function while iteratively updating the signal support using the gradient of the objective function, with the additional constraints that the signal is non-negative. The objective function may be a constrained nonlinear multivariable function. It is "sparsity-constrained" in the sense that sparser solutions have lower cost. For example, the objective function may have the following form: min _x \\ \Ax\ ² - y\\ ₂

s. t. ||x|| o < 5

supp(x) G {1, 2, ... , n}, x G R ^N

In this expression, y is a measured signal, A is a product of Gauss-Hermite and Discrete Fourier Transform matrices, and || || ₀ stands for the zero-"norm," i.e., the number of nonzero elements in the matrix.

The sparsity-constrained search performed at step 66 may be implemented in various ways. One example, using the MatLab Optimization Toolbox and Constrained Nonlinear Optimization function (fmincon) with an interior-point algorithm, is presented below in an Appendix. Alternative implementations will be apparent to those skilled in the art and are considered to be within the scope of the present invention.

After finding the optimal sparse set of basis functions and coefficients, processor 52 outputs the resulting waveform, at a pulse profile output step 68. In the present embodiment, the waveform provides a profile of the input laser pulse intensity as a function of time. Fig. 5 is a plot that schematically illustrates the intensity waveform of the pulse whose spectrum is shown in Fig. 2, following reconstruction by the method of Fig. 4, in accordance with an embodiment of the present invention. The time scale is in femtoseconds, while the intensity is in arbitrary units. This waveform includes structure that is reconstructed with a resolution less than half of the total pulse width (roughly 50 fs FWHM), and is less than one quarter of the width of the peak in the autocorrelation trace (Fig. 3). For example, the subsidiary peak occurring just prior to the main peak of the waveform has a width of about 20 fs. The waveform reconstructed and shown in Fig. 5 agreed with a FROG-based measurement of the same pulse train to within the measurement noise level. The sparsity level of the reconstruction found at step 66 and illustrated in Fig. 5, however, was eight, i.e., only eight of the functions in the Gauss-Hermite basis had non-zero coefficients. Specifically, the basis functions ip _n , _m ,q with indices { ip _{4 Q} , i/> ₅ ,i,o> ^5,3,0» 'e.i.o. ^6,3,0 > ^7,1.0 » ^7,3,0» ^7,5,0 } had corresponding coefficients {0.49, 0.78, 0.19, 0.23, 0.08, 1.87, 0.05, 0.73 }.

Fig. 6 is a plot showing cumulative distribution functions of reconstructed pulses as a function of sparsity level, using several different autocorrelation support durations, in accordance with an embodiment of the present invention. (The sparsity level, as explained above, is the minimal-number of the Gauss-Hermite functions, selected from the frame of 180 functions, required to express the intensity structure of the pulse waveform to within the target accuracy.) The plot is based on sets of simulated laser pulses having the power spectrum shown in Fig. 2 and fifth-order polynomial spectral chirps with random coefficients. Each set consisted of 5000 pulses, with a different range of support.

Traces 70, 72, 74 and 76 in Fig. 6 show the cumulative distribution functions of the sets of pulses with autocorrelation supports of 200 fs, 400 fs, 600 fs, and 800 fs, respectively. Arrows 80, 82, 84 and 86 respectively indicate the sparsity levels at which the cumulative distribution levels reach one, meaning that all of the pulses in the corresponding set could be represented by this number of Gauss-Hermite functions. For comparison, a point 78 is marked to show the sparsity level of a transform-limited pulse, having duration 42 fs and autocorrelation support of 160 fs. Fig. 6 thus demonstrates that given a reasonable choice of autocorrelation support, relative to the width of the pulse to be reconstructed, the sparsity level of the reconstructed waveform is no more than fifty, and will be no more than seventeen in nearly all cases.

As explained earlier, the results illustrated in the figures were computed by reconstructing the pulse waveforms subject to a sparsity constraint. The inventors have found that when sparsity is not constrained in the reconstruction process, the resulting pulse reconstructions have both much higher sparsity levels and much higher levels of error.

Although the embodiment described above makes use of the Gauss-Hermite functions as a convenient basis for representation of ultra-short laser pulses, other sets of functions, as are known in the art, may alternatively be used in constructing bases for this purpose. Further alternatively, depending on the actual underlying structure of the pulses in question, the methods described above may be applied, mutatis mutandis, in reconstruction of waveforms using sparse representations in terms of other functional bases. Given additional spectral information, the present sparsity-based techniques may be extended to reconstruct pulse amplitude and phase profiles, and not only the intensity waveform as described above Furthermore, although the disclosed embodiments refer to optical laser pulses, the principles of the present invention may similarly be applied in reconstructing the waveforms of pulses of other sorts.

It will thus be appreciated that the embodiments described above are cited by way of example, and that the present invention is not limited to what has been particularly shown and described hereinabove. Rather, the scope of the present invention includes both combinations and subcombinations of the various features described hereinabove, as well as variations and modifications thereof which would occur to persons skilled in the art upon reading the foregoing description and which are not disclosed in the prior art. APPENDIX - GESPAR IMPLEMENTATION

The software code presented below uses the MatLab Optimization Toolbox and Constrained Nonlinear Optimization function (fmincon), with an interior-point algorithm. The inventors have found this code to achieve efficient phase retrieval, which also leads to good recovery performance. This algorithmic method, which is termed GESPAR (GrEedy Sparse PhAse Retrieval), as described in the reference by Schechtman et al. that is cited below, is based on a fast 2-opt local search method and non-linear optimization of a sparsity-constrained formulation.

The main idea underlying GESPAR is a local-search method based on a damped Gauss- Newton method, with iterative update of the signal support using the gradient of the objective function. In the present short-pulse application, the inventors used positivity as prior information, in addition to sparsity. Accordingly, the damped Gauss-Newton step in the original GESPAR algorithm was replaced by a more general step involving minimization of a constrained nonlinear multivariable function incorporating the positivity constraint. References:

Shechtman, Y., Beck, A. & Eldar, Y.C., "GESPAR: Efficient Phase Retrieval of Sparse Signals," IEEE Transactions on Signal Processing 62, pages 928-938 (2014).

C. Papadimitriou, K. Steiglitz, Combinatorial optimization: algorithms and complexity (Dover publications, 1998).

R.H. Byrd, M.E. Hribar, J. Nocedal, "An Interior Point Algorithm for Large-Scale Nonlinear Programming," SI AM Journal on Optimization 9, pages 877-900 (1999).

MATLAB CODE LISTING

%% main function

% cd(fileparts(mfilename('main_vl.m')))

data = load('pulse_data.mat');

load('AC_data.mat');

p_ref = data.pulse;

Nt = length(p_ref);

clear data;

%% Dictionary

D = zeros(128,60);

Find = 1;

Ng = 5;

forjn = l : 10

D(:,Find:Find+Ng) = GH_mtx(15,-2.6+jn/3,Ng,Nt)';

Find = Find+Ng+1;

end

%%

F = dftmtx(Nt);

M = F*D;

c = abs(fft(AC));

%%

iterMAX = 20;

tol = 2.128e-2;

k=l;

flag=0;

while k<20 && flag == 0 [ίΜίη,χΒ] = GESPAR_M(M, c, D, k, zeros(Nt,l), iterMAX, tol);

fprintf('\n Sparsity= %d Error = %d Tol = %d\n', k, fMin,tol); if fMin<tol

flag=l;

cprintf( '*blue', '\nSolution found :) Sp=%d\n', k);

else

cprintf('err','\nNo solution found increasing sprsity :( ==> s=%d\n\n',k+l);

end

pause(2);

k = k+l;

end

%%

sol = bestMatch(D*xB,p_ref);

time = linspace(-300, 300, Nt); figure 1 = figure; % Create axes

axesl = axes('Parent',figurel,'FontSize',16);

box(axesl,'on');

hold(axesl,'air);

plot(time, p_ref/max(p_ref),'LineWidth',3,'Color',[l 0 0],'DisplayName','FROG') hold all

plot(time, sol/max(sol),'LineWidth',3,'Color',[0 0 l],'LineStyle','-','DisplayName','AC') xlabel(Time [fsec]','FontSize',16);

ylabel('Intensity [a.u.] ','FontSize', 16) ;

legend(axesl,'show')

function [fMin,xB] = GESPAR_M(M,c,Dic,k,x_ref,iterMAX,tol)

[nn, ~] = size(M);

iter = 1 ; flag = inf;

while iter<iterMAX && flag>tol

[fMin,x_n] = GreedysparseRec(M,c,k,l:nn,Dic);

% x_temp = bestMatch(Dic*x_n,x_ref);

x_temp = x_n;

if f in<flag

flag = fMin;

xB = x_temp;

end

fprintf('Iter: %d Error = %d Best Error = %d\n', iter, fMin.flag);

iter = iter+1;

end end function [fMin,x_n] = GreedysparseRec(M,c,k,measurementSet,Dic) verbose=0;

% Greedy sparseRec Finds a k sparse solution x to c=IMxl ^A 2, using only the % measurements in measurementSet

Mind = 8;

Pind = 8;

randomizeReplacement= 1 ;

[m,n]=size(M);

p=randperm(n);

supp=p(l:k);

iterations=200; % GN iterations

ind=0;

for kk=measurementSet

ind=ind+l;

% A{ind}=M(kk,:)'*M(kk,:);

A{ind}=real(M(kk,:)'*M(kk,:)) + imag(M(kk,:)'*M(kk,:));

end

M=M (measurementSet, : ) ; c=c(measurementSet) ;

% w=(l+(rand(length(measurementSet),l)<0.5))'; % Use random weights

w=ones(size(measurementSet));

% figure ;plot(real(x_k))

% [x_k,iterationEr]=GN(supp,A,c,n,randn(k,l),iterations,M,w); % Initial guess - Gauss Newton

[x_k]=GN_PS(M,supp,c,Dic);

% figure;plot(diff(-iterationEr));pause; suppTemp=supp;

fMin=WG_cost(M,c,x_k,w);

it=0;

while(l)

it=it+l;

%% Main iteration

[junk,idx]=sort(abs(x_k(supp)));

supp=supp(idx); % Sorting supp from min(abs(x_k)) to max

fGrad=WGf_QU_G_Gradient(A,c,x_k);

offSupp=setdiff(l:n,supp);

[junk,idx]=sort(-abs(fGrad(offSupp)));

offSupp=offSupp(idx);

if (randomizeReplacement)

if (length(supp)>Mind)

supp=[supp(randperm(Mind)),supp(Mind+l :end)] ;

else

supp=supp(randperm(length(supp)));

end

if (length(offSupp)>Pind)

offSupp=[offSupp(randperm(Pind)),offSupp(Pind+l :end)] ;

else

off S upp=off S upp(randperm(length(of f S upp)) ) ;

end

end pSupp=l:length(supp);

pOff S upp= 1 : length(of f S upp) ;

improved=0;

for ilnd=l:lengfh(supp)

i=supp(pSupp(iInd)); %Index to remove

for jlnd=l :min(Pind,length(pOffSupp))

j=offSupp(pOffSupp jInd)); % Index to insert

%% Check replacement

suppTemp=supp ;

suppTemp((suppTemp==i))=j ;

%% Solve GN with given support

% xTemp=GN(suppTemp,A,c,n,x_k(suppTemp),iterations,M,w);

[xTemp] =GN_PS (M,suppTemp,c,Dic);

fTemp=WG_cost(M,c,xTemp,w);

if fTemp<fMin

if (verbose)

fprintf('it: %d Replaced d with %d f= %3.3f\n',it, i,j,fTemp); suppTemp

end

x_k=xTemp;

x_n=x_k;

supp=suppTemp;

improved=l;

fMin=fTemp;

if fTemp<le-3

if (verbose)

fprintf('* * ****************************************5 ucces s ! , iteration=%d\n',it) ;end return;

end break;

end

end if (improved)

break;

end

end

if (-improved)

if (verbose)

fprintf('no possible improvement - trying new initial guess\n'); end

x_n=x_k;

return

end

end

x_n=x_k; end function [x_t, fval]=GN_PS(M,supp,c,Dic)

b = zeros(size(Dic,l),l)+0.00;

x_t = zeros(size(M,2),l);

xO = rand(size(supp))';

M = M(:,supp);

f = @(x)Afun(x,M,c);

options = optimset('Display','off ,'AlgorithmVactive-set');

% [x,fval] = fminunc(f, xO, Die, b, options);

[x,fval] = fmincon(f, xO, -Dic(:,supp), b,[], [], [], [], [], options); x_t(supp) = x;

end function y = Afun(x,M,c)

y = sqrt( sum( abs(abs(M*x). ^A 2-c). ^A 2 ) );

end function out=WG_cost(A,c,x,w)

% Returns the squared 2 norm of the consistency error % out=norm(abs(A*x). ^A 2-c);

out=sqrt(sum(w'.*((abs(A*x) ^A 2-c) ^A 2)));

end function out=WGf_QU_G_Gradient(A,c,x)

% Returns gradient of cost function around current estimate x out=0;

m=length(A);

for i=l:m

out=out+4*(x'*A{i}*x-c(i))*A{i}*x;

end

end function [GH] = GH_mtx(T,tau,Ng,Nt)

% %return HG frame

time = linspace(-T,T,Nt)-tau;

p=hermipol(Ng);

for i=l:length(p)

temp = p(i,l:i);

HG(i,:) = polyval(temp,time).*exp(-(time. ^A 2)/2);

HG(i,:) = HG(i,:)/sqrt(HG(i,:)*HG(i,:)');

end

GH = HG; function p=hermipol(n) k=2;

P(U)=1 ;

p(2,l :2)=[2 0];

for k=2:n

p(k+ 1 , 1 :k+ 1 )=2* [p(k, 1 :k) 0] -2* (k- 1 )* [0 0 p(k- 1 , 1 :k- 1 )] ; end