classdef FarrowDesignUtils
    % Class containing functions for the LTE sample rate conversion example.
    
    methods (Static)
        
        function w = designFilterPrototype(farrow)
            % Expects a structure farrow with the following elements:
            %      FsIn: Input sampling rate.
            %     FsOut: Output sampling rate.
            %      Fsig: Bandwidth of each sideband of the signal of interest.
            %   numTaps: Filter length as a multiple of input sample times.
            %       sps: Oversampling factor of prototype filter.
            
            if farrow.FsOut < farrow.FsIn
                % If FsOut is less than FsIn, adjust the filter to prevent aliasing
                % of the LTE signal spectrum at the new sampling frequency
                Fdelta = farrow.FsIn-farrow.FsOut;
            else
                % Otherwise, the filter
                Fdelta = 0;
            end
            
            % To create something close to a Nyquist filter, Fpass and Fstop must be related by the
            % equations below. To prevent aliasing at the new sampling rate,
            % increase Fpass by Fdelta, which has the effect of decreasing
            % Fstop by Fdelta.
            
            % Max LTE bandwidth is 20 MHz; +/-10MHz on each side of the carrier.
            Fpass  = farrow.Fsig + Fdelta;
            F6db   = farrow.FsIn/2;             % Aim for 6 dB midway through
            Fstop  = farrow.FsIn-Fpass;
            Fsover = farrow.FsIn * farrow.sps;  % Sampling rate of oversampled filter
            order  = farrow.sps*farrow.numTaps; % Order of prototype filter
            
            Rp  = (10^(0.01/20)-1)/2;  % 0.01 dB peak-to-peak ripple
            Rst = 1e-3;                % 60 dB stopband attenuation
            
            w = firgr(order,...
                [0,Fpass,F6db,F6db,Fstop,Fsover/2]/(Fsover/2),...
                [1 1 0.5 0.5 0 0],...
                [Rp Rp Rst]);
            
            % Multiply coefficients by oversampling factor to compensate for gain.
            w = farrow.sps * w;
            
        end
        
        
        function A = generatePolynomialCoefficients(farrow)
            
            % Compute polynomial coefficients for piecewise cubic approximation of the
            % prototype filter. Each slice spans sps (samples per slice) samples.
            
            order = 3;
            proto = farrow.prototype;
            sps = farrow.sps;
            
            % Calculate number of slices, L.
            L = floor(length(proto)/sps);
            r = rem(length(proto),L);
            
            if (r ~= 1) || (rem(L,2)~=0)
                error(['Protoype filter length must be equal to L*ovsf+1 where L is an even integer.' newline ...
                    'Prototype filter length is ' num2str(L) '*ovsf+' num2str(r) '.']);
            end
            
            A = zeros(order+1,L);
            
            % Fit L curves to raised cosine
            for k = 0:L-1
                x        = ((0:sps)./sps);
                y        = proto(sps*k+1:sps*k+1+sps);
                A(:,k+1) = flipud(polyfit(x,y,order).');
            end
            
        end
        
        function [s,t] = evaluateApproximation(polynomials,sps)
            
            numTaps = size(polynomials,2);
            
            s = zeros(numTaps*sps+1,1);
            t = zeros(numTaps*sps+1,1);
                        
            for k=1:numTaps
                x = (0:sps-1)/sps;
                s((k-1)*sps + (1:sps)) = polyval(flipud(polynomials(:,k)),x);
                t((k-1)*sps + (1:sps)) = x + (k - 1 - numTaps/2);
            end
            
            % last point
            s(end) = polyval(flipud(polynomials(:,numTaps)),1);
            t(end) = numTaps/2;
            
        end
        
        function [y, t] = convertSampleRate(farrow,x)
            
            A = farrow.polynomials;            
            R = farrow.FsOut/farrow.FsIn;
            stride = 1/R;
            
            numFIRs = size(A,1);
            r       = zeros(length(x),numFIRs);
            
            outputLength = round(length(x)*R);
            outputTimes  = 0:stride:outputLength*stride-1;
            
            % Compute FIR outputs
            for m=1:numFIRs
                r(:,m) = filter(A(m,:),1,x);
            end
            
            % Compute the mult-add chain output for each output sample
            y = zeros(outputLength,1);
            
            for n=1:outputLength
                kz     = floor(outputTimes(n));
                rho    = outputTimes(n)-kz;
                rhoVec = [1 rho rho^2 rho^3];
                y(n)   = rhoVec * r(kz+1,:).';
            end
            
            t = outputTimes;
            
        end
        
        %---------------------------------------------------------------------------------------------------
        % Helper for plotting an LTE signal region image
        %---------------------------------------------------------------------------------------------------
        
        function plotSignalImages(FsOut)
            
            ax        = axis;
            ymin      = ax(3) - 10;
            xtemplate = [0 10e6 10e6 FsOut-10e6 FsOut-10e6 FsOut]/1e6;
            ytemplate = [-3 -3 ymin ymin -3 -3];
            
            nmax = floor(max(ax(2))/(FsOut/1e6));
            
            % Image centre frequency indices
            fcs = FsOut * (-nmax:nmax-1) / 1e6;
            
            x = repmat(xtemplate,1,length(fcs));
            y = repmat(ytemplate,1,length(fcs));
            
            % Add centre frequency offsets to x
            for k=0:length(fcs)-1
                x((1:6)+6*k) = x((1:6)+6*k) + fcs(k+1);
            end
            
            plot(x,y,'r');
            
            % Restore original axis
            axis(ax);
            
        end
        
        
    end
    
end

