% Polyphase FIR with symmetric impulse response
% Markus Nentwig, 12.10.2012
%
% A FIR filter can be implemented using one multiplication per input sample
% and independent coefficient.
% This allows to exploit a symmetric impulse response.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

function pphase()
    close all;
    
    % **************************************************************
    % A symmetric impulse response
    % **************************************************************
    h = [1 3 5 7 11 13 11 7 5 3 1];    

    % one input sample generates nUps output samples
    nUps = 3;

    % all coefficients behind the center (11) are duplicates
    nIndepCoeff = (numel(h)+1)/2; % which is 6

    % extract the independent coefficients
    hTrunc = h(1:nIndepCoeff); % which is [1 3 5 7 9 11]

    % Over how many input samples does the impulse response extend?
    memDepth = ceil(numel(h) / nUps); % which is 4

    % test signal length    
    nSamplesIn = 50;
    nSamplesOut = nSamplesIn * nUps;
    
    if false
        % generate impulse test signal
        sig = zeros(1, nSamplesIn);
        sig(10) = 1;
    else
        % generate arbitrary test signal
        sig = rand(1, nSamplesIn);
    end

    % **************************************************************
    % reference: convolution with impulse response
    % **************************************************************    

    % insert (nUps-1) zeros for each input sample
    sig1 = [sig; zeros(nUps-1, nSamplesIn)];
    sig1 = sig1(:) .';
    
    % convolve with impulse response
    sig1 = conv(sig1, h);
    sig1 = sig1(1:nSamplesOut);
    
    % **************************************************************
    % conventional polyphase implementation
    % **************************************************************    
    sig2a = [];
    x = zeros(1, memDepth);
    for ix = 1 : nSamplesIn
        % shift input sample into delay line
        x = [sig(ix) x(1:end-1)];

        % each input sample generates three output samples
        % 10 multiplications / input sample
        % 8 additions / input sample
        y1 = x(1) * h(1) + x(2) * h(4) + x(3) * h(5) + x(4) * h(2);
        y2 = x(1) * h(2) + x(2) * h(5) + x(3) * h(4) + x(4) * h(1);
        y3 = x(1) * h(3) + x(2) * h(6) + x(3) * h(3); % -> (x(1)+x(3))h(3)+...
        
        sig2a = [sig2a, y1, y2, y3];
    end
    
    % **************************************************************
    % efficient implementation stores and reuses sample-coefficient 
    % products
    % **************************************************************    
    % symbols from blog article:
    % A = mem(1, 1)
    % B = mem(2, 4)
    % C = mem(3, 5) etc 
    % first (row) index: delay
    % second (column) index: coefficient
    mem = zeros(memDepth, nIndepCoeff);
    sig2b = [];
    for ix = 1 : nSamplesIn
        x = sig(ix);
        % 6 multiplications/input sample (size of hTrunc)       
        prodRow = hTrunc * x; 

        % shift each product into its own delay line
        % 6 columns
        % 8 additions / input sample
        mem = [prodRow; ...
               mem(1:end-1, :)];
        y1 = mem(1, 1) + mem(2, 4) + mem(3, 5) + mem(4, 2);
        y2 = mem(1, 2) + mem(2, 5) + mem(3, 4) + mem(4, 1);
        y3 = mem(1, 3) + mem(2, 6) + mem(3, 3);
        
        sig2b = [sig2b, y1, y2, y3];
    end 
    
    % **************************************************************
    % plot outputs, should be identical
    % **************************************************************    
    figure(); hold on;
    h = plot(1:nSamplesOut, sig1, 'kx-'); set(h, 'lineWidth', 3);
    h = plot(1:nSamplesOut, sig2a, 'b+-'); set(h, 'lineWidth', 2);
    h = plot(1:nSamplesOut, sig2b, 'g+-'); set(h, 'lineWidth', 1);
    legend('convolution', 'conventional implementation', 'polyphased implementation');
end

