function [RSP,ISP]=SRf(SP,SS)
%RECOVERY PART
w0=SP(1);
%Make LF extrapolation
DSP=diff(SS)./diff(SP);
rs=DSP(1);%must be negative
r=SS(1);
alp=rs/r/w0-log(r)/(w0*w0);
bet=2*log(r)/w0-rs/r;
x=1:(w0-1);
LFpart=exp(alp.*x.*x+bet.*x-alp-bet);
%Make interpolation
ISP = pchip(SP,SS,w0:max(SP));
%Make HF extrapolation

w=1:(2^15-length(ISP)-length(LFpart));
wmax=max(SP);
rosh=DSP(end);

  alp=(rosh/ISP(end))*wmax;
        c=0;
        while (alp>0)
            c=c+1;
          rosh=  (SS(end)-SS(end-c))/(SP(end)-SP(end-c));
          alp=(rosh/ISP(end))*wmax;
        end
HFpart=ISP(end)* (1+w/wmax).^(alp);
HFpart(HFpart<eps)=eps;

%Connect all together
ISP=[LFpart,ISP,HFpart];
%Mirror image of the spectrum
RSP=[ISP,fliplr(ISP)];