function plant = freblkData2Frd(data,w,units)
%% Offline FFT estimation based on the logged data from FRE block
%
%   data: a structure of 4 fields: d/u/y/ready, each a timeseries object.
%   All of them must be uniformed sampled with the same Ts.
%
%   w and units: same frequencies and units used in the FRE block.
%
%   usefiltering is used by sine stream FFT.

%   Author: R. Chen
%   Copyright 2018 The MathWorks, Inc.

%% process units
usefiltering = true;
w = unitconv(w,units,'rad/s');
nw = length(w);
%% get time series
dts = data.Perturbation;
uts = data.PlantInput;
yts = data.PlantOutput;
rts = data.Ready;
Ts = yts.TimeInfo.Increment; % use yts as the truth because it is always sampled properly
dupidx = find(diff(rts.Time)<Ts/2);
rts = delsample(rts,'Index',dupidx);
dupidx = find(diff(dts.Time)<Ts/2);
dts = delsample(dts,'Index',dupidx);
dupidx = find(diff(uts.Time)<Ts/2);
uts = delsample(uts,'Index',dupidx);
%% determine which mode is used from the ready signal (# of rising edges)
IsNotSinestreamMode = (nw>1) && sum(diff(rts.Data)>0)==1; % only has one edge
%% Estimate FRD from d, u and y signals
if IsNotSinestreamMode
    %% get measured signals (d start from 0)
    idx = rts.Data==1;
    d = dts.Data(idx);
    u = uts.Data(idx);
    y = yts.Data(idx);
    %% get FFT
    NFFT = sum(idx);
    if rem(NFFT,2)
        % Odd
        freq = 2*(pi/Ts)*linspace(0,1,NFFT+1);
        last_in = (NFFT+1)/2;
        freq = freq(1:last_in);
    else
        % Even
        freq = (pi/Ts)*linspace(0,1,NFFT/2+1);
    end

    %% create filter window
    pd = 2*pi/(NFFT-1);
    filterWindow = [0.5*ones(NFFT-1,1)-0.5*cos(pd*(0:NFFT-2)'); 0];
    %% apply filter window to signal
    fd = fft(d.*filterWindow,NFFT);
    fu = fft(u.*filterWindow,NFFT);
    fy = fft(y.*filterWindow,NFFT);
    d2u = fu./fd;
    d2y = fy./fd;
    u2y = d2y./d2u;
    %% generate FRD from interpolation
    frdFFT = frd(u2y(1:length(freq)), freq, Ts);
    % Only resample the FRD object for results obtained using PRBS input
    % signals with NumPeriods>1
    if isfield(data,'Info') && isfield(data.Info,'PRBSNumPeriods')
        idxInfo = idx(end);
        NumPeriods = data.Info.PRBSNumPeriods.Data(idxInfo);
        if NumPeriods > 1
            ResponseData = frdFFT.ResponseData(:,:,1:NumPeriods:end);
            Frequency = frdFFT.Frequency(1:NumPeriods:end);
            frdFFT = frd(ResponseData, Frequency, Ts);
        end 
    end
    plant = interp(frdFFT, w);
    plant = chgFreqUnit(plant,units);
else
    %% find start, end and samples for each frequency
    startEstimationPoints = find(diff(rts.Data)>0);
    endEstimationPoints = find(diff(rts.Data)<0);
    % if last experiment not completed, end point is the last measurement
    if length(endEstimationPoints)<length(startEstimationPoints)
        endEstimationPoints = [endEstimationPoints;length(rts.Data)];
    end
    beginFreqPoints = [0;endEstimationPoints(1:end-1)] + 1;
    sampsPerPeriod = round(2*pi/Ts./w);
    %% compute frequency response
    u2y = zeros(nw,1); 
    selectedW = zeros(nw,1);
    for ct = 1:nw
        if usefiltering
            % design a bandpass filter
            filt = frest.Sinestream.designFIRFilter(Ts,sampsPerPeriod(ct));
            % extract the full frequency data including settling periods
            uThisFreq = uts.Data(beginFreqPoints(ct):endEstimationPoints(ct));
            yThisFreq = yts.Data(beginFreqPoints(ct):endEstimationPoints(ct));
            dThisFreq = dts.Data(beginFreqPoints(ct):endEstimationPoints(ct));  
            % filter
            uThisFreq = filter(filt,1,uThisFreq);
            yThisFreq = filter(filt,1,yThisFreq);
            dThisFreq = filter(filt,1,dThisFreq);
            % cut out the transient data
            uThisFreq = uThisFreq((startEstimationPoints(ct)-beginFreqPoints(ct)+1):end);
            yThisFreq = yThisFreq((startEstimationPoints(ct)-beginFreqPoints(ct)+1):end);
            dThisFreq = dThisFreq((startEstimationPoints(ct)-beginFreqPoints(ct)+1):end);
        else
            % cut out the transient data
            uThisFreq = uts.Data(startEstimationPoints(ct):endEstimationPoints(ct));
            yThisFreq = yts.Data(startEstimationPoints(ct):endEstimationPoints(ct));
            dThisFreq = dts.Data(startEstimationPoints(ct):endEstimationPoints(ct));  
        end
        % compute FFT
        NFFTThisFreq = numel(uThisFreq);    
        dfft = fft(dThisFreq(:),NFFTThisFreq);
        d2uThisFreq = fft(uThisFreq(:),NFFTThisFreq)./dfft;
        d2yThisFreq = fft(yThisFreq(:),NFFTThisFreq)./dfft;
        u2yThisFreq = d2yThisFreq./d2uThisFreq;
        %% generate FRD from interpolation
        if rem(NFFTThisFreq,2)
            % Odd
            freq = 2*(pi/Ts)*linspace(0,1,NFFTThisFreq+1);
            last_in = (NFFTThisFreq+1)/2;
            freq = freq(1:last_in);
        else
            % Even
            freq = (pi/Ts)*linspace(0,1,NFFTThisFreq/2+1);
        end
        frdFFT = frd(u2yThisFreq(1:length(freq)), freq, Ts);
        plant = interp(frdFFT, w(ct));
        u2y(ct) = plant.ResponseData;
    end
    %% generate frd plant
    plant = frd(u2y,unitconv(w,'rad/s',units),'Ts',Ts,'FrequencyUnit',units);
end


