diff --git a/SDP/SPP/MATLAB/polyphase_filter/partition.m b/SDP/SPP/MATLAB/polyphase_filter/partition.m
new file mode 100755
index 0000000000000000000000000000000000000000..4cd4bd4d9fda4f098ea29baf4d8d130cd09f1267
--- /dev/null
+++ b/SDP/SPP/MATLAB/polyphase_filter/partition.m
@@ -0,0 +1,54 @@
+function [indices] = partition(n,m,l,option)
+%
+% function [indices] = partition(n,m,l,option)
+%
+% input
+% n - length of the signal
+% m - length of a block
+% l - length of overlapping section (default=0)
+% option - 'set' [default] returns the entire set of indices
+% 'range' returns a list of ranges in pairs of [first,last]
+% output:
+%
+% In case option='set' the returned value indices is a floor( (n-m) /(m-l))+1 by m array
+% each row contains the indices for one block of data, the rows are consequetive, possibly
+%` overlapping windows, e.g. partition(8,4,2) returns:
+%
+% 1 2 3 4
+% 3 4 5 6
+% 5 6 7 8
+%
+% In cse option='range' the returned value indices is a floor( (n-m) /(m-l))+1 by 2 matrix
+% with indices(:,1) are the first indices of the ranges and indices(:,2) are the last indices
+% of the ranges, e.g. partition(8,4,2,'range') returns:
+%
+% 1 4
+% 3 6
+% 5 8
+%
+% See also: reordering,
+%
+% (C) 2002 Astron, by M.van Veelen
+
+if nargin < 4
+ option='set';
+end ;
+
+if nargin<3
+ l=0;
+end;
+
+
+n_indices = floor( (n-m) /(m-l))+1;
+
+if strcmp(option,'range')
+ indices=zeros(n_indices,2) ;
+ indices(:,1) = [1:(m-l):n-m+1 ]';
+ indices(:,2) = indices(:,1)+m-1;
+else
+ indices=zeros(n_indices, floor(m) ) ;
+ for i=0:n_indices-1
+ indices(i+1,:)=i*(m-l)+[1:m] ;
+ end;
+end;
+
diff --git a/SDP/SPP/MATLAB/polyphase_filter/polyphase.m b/SDP/SPP/MATLAB/polyphase_filter/polyphase.m
new file mode 100755
index 0000000000000000000000000000000000000000..a9684857361be3cb5abaee3ee9177cafbcccd035
--- /dev/null
+++ b/SDP/SPP/MATLAB/polyphase_filter/polyphase.m
@@ -0,0 +1,59 @@
+function [Yall,Yall_test] = polyphase(x,K,L,M,option)
+%
+% function [Y] = polyphase(x,K,L,M,option)
+%
+% K - the number of channels
+% L - filter length
+% M - the downsampling rate
+% option - 'normal' [default]
+% 'shift' applies fftshift on the output of the filterbank
+% 'centre' centers the signals and recombines the output accordingly
+%
+% Compute the polyphase filter coefficients p and applies them the the data x, i.e. the
+% data is decimated by a factor M and divided into K channels. The prototype filter can
+% be supplied through h, by default fircls1 is used to derive the prototype filter window.
+% The coefficients care derived from this prototype windows and convolved with x, hence
+% a number of output vectors are computed which can be transformed into the real filter
+% outputs using an DFT. The polyphase filter structure comprises of two stages:
+%
+% STAGE 1: Convolution
+% STAGE 2: Fourier Transform
+%
+% The both stages are performed. The output is an N X K matrix with a % channel in each columm,
+% providing K columns; the length N, is floor((length(x))/K)-L+1, in case the data is critically
+% sampled (K=M). Theoretical background can be found in
+%
+% [1] Multirate Digital Signal Processing. Ronald. E. Crochiere & Lawrence R.Rabiner,1983.
+% ISBM 0-13-605162-6, Prentice Hall Inc., Englewood Cliffs, New Jersey 07632.
+%
+% See also: reordering, partition, prefilter, polyphase_coefficients
+%
+% (C) 2002 Astron, by M.van Veelen
+%
+
+if nargin < 5 ; option='normal' ; end ;
+
+% Settings for finite word-length coding
+
+% INITIALIZE: compute the filter prototype coefficients of the
+centre='no' ; if strcmp(option,'centre') ; centre='yes' ; end ;
+
+% STAGE 1: Pre-filter by convolution with polyphase filter coefficients
+ym=prefilter(x,K,L,M,centre);
+
+% STAGE 2: DFT the output of the pre-filter to obtain the channels
+Yall=1/K * fft(ym)' ;
+
+if strcmp(option,'shift') | strcmp(option,'centre') ;
+ Yall= [ Yall(:,floor(K/2):-1:1) , Yall(:,K:-1:floor(K/2)+1) ];
+end ;
+
+% Normalization of the input signals in order to replace the signal in the range [-1,+1]
+% x = x / max(x);
+% x = double(uencode(x, quant_signal, 1, 'signed'));
+% x = x / max(x);
+
+
+
+
+
diff --git a/SDP/SPP/MATLAB/polyphase_filter/polyphase_coefficients.m b/SDP/SPP/MATLAB/polyphase_filter/polyphase_coefficients.m
new file mode 100755
index 0000000000000000000000000000000000000000..9deec2b3aae263af739a911390779185c9ed9717
--- /dev/null
+++ b/SDP/SPP/MATLAB/polyphase_filter/polyphase_coefficients.m
@@ -0,0 +1,112 @@
+function [p,h,p_test] = polyphase_coefficients(K,L,M,option,centre,h)
+%
+% function [p,h] = polyphase_coefficients(K,L,M)
+%
+% input:
+% K - the number of channels
+% L - filter length
+% M - the downsampling rate (default M=K)
+% option - 'coefficients' returns the coefficients of the polyphase (default)
+% 'indices' returns the indices into the prototype filter window
+% centre - if centre='yes' the output is shifted to the centre (default centre='no')
+% h - replaces the default fircls1 filter window with the one specified
+%
+% output:
+% p - the coefficients of the polyphase in case option = 'coefficients'
+% indices into the prototype filter window in case option = 'indices'
+% h - the prototype filter window
+%
+% The filter coefficients can be applied directly to the data, however it is not necessary
+% to compute the filterbank coefficients p, as the following example will show. The filter
+% coefficients of the polyphase structure can be optained in two ways (e.g. take K=16,L=8,M=16):
+%
+% [p,h] = polyphase_coefficents(K,L,M) ;
+% [i,h] = polyphase_coefficient(K,L,M,'indices') ;
+% h(i) == p
+%
+% The use of indices computation allows you to easily replace the prototype filter window.
+% The current implementation uses firls, in the following example we consider fir1:
+%
+% [i,h] = polyphase_coefficient(K,L,M,'indices') ;
+% alpha = 1 ;
+% h=fir1(K-1,L/K * alpha) ;
+% h = resample(h, K * L, K);
+% h = h / sqrt(h * h'); % nomalisation
+% h = h / max(h);
+% p=h(i) ;
+%
+% The same result can be achieved by directly providing a prototype filter window, e.g.
+%
+% hfir=fir1(K-1,L/K * alpha) ;
+% p=polyphase_coefficients(16,12,16,'coefficients','no',hfir) ;
+% [i,h]=polyphase_coefficients(16,12,16,'coefficients','no',hfir) ;
+%
+% Note that the returned filter windows is resampled so h != hfir in the example above,
+% and hfir(i) is NOT a valid indexing of hfir, but h(i)==p.
+%
+% See also: reordering, partition, prefilter, polyphase, firls
+%
+% (C) 2002 Astron, by M.van Veelen, Version 1.2
+
+p=zeros(K,L) ;
+
+if nargin < 3 ; M=K ; end ;
+if nargin < 4 ; option = 'coefficients' ; end ;
+if nargin < 5 ; centre = 'no' ; end ;
+
+if nargout > 2
+ p_test=zeros(K,L) ;
+end ;
+
+% Quantisation settings (if any)
+h_nq = 0 ; % Number of quantisation level for the prototype filter windows h
+
+% Compute the prototype filter window
+if nargin < 6
+ alpha =1 ; % alpha = M/(K);
+ h = fircls1(K-1, L / K * alpha, 0.01, 0.001);
+end;
+
+ h = resample(h, length(h) * L, K) ; % version 1.0
+ % h = resample(h, length(h) * floor(L * K/M), K) ; % version 1.1
+
+ if strcmp(centre, 'yes') ;
+ h=exp(2*sqrt(-1)*pi*L/2*(1:1:length(h))/length(h)).*(K * L * h);
+ else
+ % tryout for the first phase
+ h=exp(-2*sqrt(-1)*pi*L/2*(1:1:length(h))/length(h)).*(K * L * h);
+ end;
+
+ h = h / sqrt(h * h'); % nomalisation
+ h = h / max(h);
+
+if h_nq > 0
+ % Quantization and normalization of the filter in order to avoid a bias
+ h = double(uencode(h, quant_h, 1, 'signed'));
+ h = h / max(h);
+end;
+
+
+% p=partition(K*L,K)'; % version 1.0
+% p=partition(K*floor(L*K/M),K)' ; % version 1.1
+p=partition(K*L,K,K-M)'; % version 1.2
+
+if ~strcmp(option, 'indices')
+ p = h(p) ;
+end
+
+if nargout > 2
+ for n = 1 : K
+ for i = 1 : L
+ p_test(n, i) = h((i - 1) * K + n);
+ end ;
+ end
+end
+
+%fh = fopen('c:\temp\h.dat', 'w+');
+%fwrite(fh, SubFilter, 'double');
+%fclose(fh);
+
+
+
+
diff --git a/SDP/SPP/MATLAB/polyphase_filter/polyphase_compare.m b/SDP/SPP/MATLAB/polyphase_filter/polyphase_compare.m
new file mode 100755
index 0000000000000000000000000000000000000000..42cf0b44ee957f9f3893bf58f20a729c5636cde1
--- /dev/null
+++ b/SDP/SPP/MATLAB/polyphase_filter/polyphase_compare.m
@@ -0,0 +1,136 @@
+function [ctk,ctd,ra,rm] = polyphase_compare(option,n_tests,scale);
+%
+% polyphase_compare(option)
+%
+% option : 'simple' just compare one run with the filterbank
+% 'channels' compare the implementations by varying the number of channnels
+% 'lenght' compare the implementations by varying the number of datapoints
+% 'full' run all the tests
+% 'snap' technicians only ;-)
+% 'snapall' technicians only ;-)
+%
+% n_tests : number of tests to run.
+%
+% scale : 'normal'
+% 'loglog' (default)
+%
+% Don't be too eager with the n_tests ... it's 2^(n_tests+offset) datapoints/channels!
+% Start with a value like 4, 8 should be computable within reasonable time on 1GHz CPUs.
+%
+% See also: reordering, partition, prefilter, polyphase_coefficients, polyphase
+%
+% (C) 2002 Astron, by M.van Veelen
+%
+
+format compact;
+if nargin < 1 ; option='simple' ; end ;
+if nargin < 2 ; n_tests = 4 ; end ;
+if nargin < 3 ; scale='loglog' ; end ;
+
+n_figures = 0 ;
+
+ra=[]; rm=[];
+
+nms=2^14;
+% x=chirp(t,50,1,150,'q');
+ dt=0.001;
+ t=0:dt:dt*nms;
+ x=chirp(t,5,1,10,'q');
+% x=sin(2*pi*t*25) + sin(2*pi*t*50);
+
+
+K=16; L=8 ;
+if strcmp(option,'simple') | strcmp(option,'snap') | strcmp(option,'full') | strcmp(option,'snapall')
+ tic;
+ p=alex_polyphase_init(L,K);
+ ra=alex_polyphase(x,16,32,32,p,K,L, floor(length(x)/L) );
+ t=toc; ctk=t;
+ figure(1); subplot(321) ;
+ imagesc(sqrt(real(ra).^2+imag(ra).^2)') ; colormap gray ; xlabel('time') ; ylabel('channel');
+ title(['Previous implementation. Computing time ',num2str(t),' [s]'] ) ;
+
+ tic
+ rm=polyphase(x,K,L,K) ; % M is set to K
+ t=toc; ctd=t ;
+ n_figures=n_figures+1 ; figure(n_figures); subplot(322) ;
+ imagesc(sqrt(real(rm).^2+imag(rm).^2)') ; colormap gray ; xlabel('time') ; ylabel('channel');
+ title(['Previous implementation. Computing time ',num2str(t),' [s]'] ) ;
+
+ if strcmp(option,'snap');
+ n_figures=n_figures+1 ; figure(n_figures);
+ imagesc(sqrt(real(rm).^2+imag(rm).^2)') ; colormap gray ; xlabel('time') ; ylabel('channel');
+ title(['Previous implementation. Computing time ',num2str(ctk),' [s]'] ) ;
+ n_figures=n_figures+1 ; figure(n_figures);
+ imagesc(sqrt(real(rm).^2+imag(rm).^2)') ; colormap gray ; xlabel('time') ; ylabel('channel');
+ title(['Previous implementation. Computing time ',num2str(ctd),' [s]'] ) ;
+ n_figures=n_figures+1 ; figure(n_figures);
+ mesh(sqrt(real(rm-ra).^2+imag(rm-ra).^2)') ; colormap gray ; xlabel('time') ; ylabel('channel');
+ title(['Difference. Computing time ',num2str(ctk-ctd),' [s]'] ) ; shading flat;
+ end ;
+end ;
+
+if strcmp(option,'full') | strcmp(option,'channels') | strcmp(option,'snapall') ;
+
+ display('running test one ...') ;
+ offset=3 ;
+ ctk=zeros(n_tests,2);
+ for n=1:n_tests;
+ K=2^(n+offset-1);
+ L=K/2;
+ tic; p=alex_polyphase_init(L,K); ra=alex_polyphase(x,16,32,32,p,K,L,1E10); ctk(n,1)=toc ;
+ tic; rm=polyphase(x,K,L,K) ; ctk(n,2)=toc ;
+ clear('p') ; clear('rm') ;clear('ra') ;
+ end;
+
+ if strcmp(scale,'loglog') | strcmp(option,'snapall');
+ if strcmp(option,'snapall') ; n_figures=n_figures+1 ; figure(n_figures); subplot(111) ;else ; subplot(313) ; end ;
+ loglog( 2.^(offset+[1:length(ctk)]), ctk(:,1),'-', 2.^(offset+[1:length(ctk)]), ctk(:,2),'--' ) ;
+ legend('Previous version', 'New version',2) ;
+ ylabel('time [s]'); xlabel('Number of channels')
+ title (['Computing time as function of number of channels (',int2str(length(x)),' data points) L=K/2' ]) ;
+ axis([ 2^(offset-1),2^(offset+n_tests+1), min(min(ctk))*0.9, max(max(ctk))*1.1 ]) ;
+ if strcmp(option,'snapall') ; saveas(gcf,['performance_channels_logscale'],'bmp'); end ;
+ end ;
+ if strcmp(scale,'normal') | strcmp(option,'snapall');
+ if strcmp(option,'snapall') ; n_figures=n_figures+1 ; figure(n_figures); subplot(111) ;else ; subplot(313) ; end ;
+ plot( 2.^(offset+[1:length(ctk)]), ctk(:,1),'-', 2.^(offset+[1:length(ctk)]), ctk(:,2),'--' ) ;
+ legend('Previous version', 'New version',2) ;
+ ylabel('time [s]'); xlabel('Number of channels')
+ title (['Computing time as function of number of channels (',int2str(length(x)),' data points) L=K/2' ]) ;
+ axis([ 2^(offset)*0.9,2^(offset+n_tests)*1.1, min(min(ctk))*0.9, max(max(ctk))*1.1 ]) ;
+ if strcmp(option,'snapall') ; saveas(gcf,['performance_channels'],'bmp'); end ;
+ end;
+
+end ;
+
+if strcmp(option,'length') | strcmp(option,'full') | strcmp(option,'snapall')
+ display('running test two ...') ;
+ ctd=zeros(n_tests,2);
+ dt=0.001; K=64; L=32 ;
+ offset=log2(K*L);
+ for n=1:n_tests;
+ nms=2^(n+offset-1);
+ t=0:dt:dt*nms;
+ x=chirp(t,5,1,10,'q');
+ tic; p=alex_polyphase_init(L,K); ra=alex_polyphase(x,16,32,32,p,K,L,1E10); ctd(n,1)=toc ;
+ tic; rm=polyphase(x,K,L,K) ; ctd(n,2)=toc ;
+ clear('p') ; clear('rm') ;clear('ra') ; clear('x') ;
+ end;
+
+ if strcmp(scale,'loglog') | strcmp(option,'snapall');
+ if strcmp(option,'snapall') ; n_figures=n_figures+1 ; figure(n_figures); subplot(111) ;else ; subplot(313) ; end ;
+ loglog( 2.^(offset+[1:length(ctd)]), ctd(:,1),'-', 2.^(offset+[1:length(ctd)]), ctd(:,2),'--' ) ;
+ legend('Previous version', 'New version',2) ; ylabel('time [s]'); xlabel('Number of samples')
+ title (['Computing time as function of number of samples (K=',int2str(K),';L=',int2str(L), ')']) ;
+ axis([ 2^(offset-1),2^(offset+n_tests+1), min(min(ctd))*0.9, max(max(ctd))*1.1 ]) ;
+ if strcmp(option,'snapall') ; saveas(gcf,['performance_datasize_logscale'],'bmp'); end ;
+ end ;
+ if strcmp(scale,'normal') | strcmp(option,'snapall');
+ if strcmp(option,'snapall') ; n_figures=n_figures+1 ; figure(n_figures); subplot(111) ;else ; subplot(313) ; end ;
+ plot( 2.^(offset+[1:length(ctd)]), ctd(:,1),'-', 2.^(offset+[1:length(ctd)]), ctd(:,2),'--' ) ;
+ legend('Previous version', 'New version',2) ; ylabel('time [s]'); xlabel('Number of samples')
+ title (['Computing time as function of number of samples (K=',int2str(K),';L=',int2str(L), ')']) ;
+ axis([ 2^(offset)*0.9,2^(offset+n_tests)*1.1, min(min(ctd))*0.9, max(max(ctd))*1.1 ]) ;
+ if strcmp(option,'snapall') ; saveas(gcf,['performance_datasize'],'bmp'); end ;
+ end;
+end ;
diff --git a/SDP/SPP/MATLAB/polyphase_filter/polyphase_demo.m b/SDP/SPP/MATLAB/polyphase_filter/polyphase_demo.m
new file mode 100755
index 0000000000000000000000000000000000000000..79b0970e0e30a5dc7b592f8b2c295a32c99bcb37
--- /dev/null
+++ b/SDP/SPP/MATLAB/polyphase_filter/polyphase_demo.m
@@ -0,0 +1,67 @@
+function polyphase_demo(option,line_style)
+%
+% polyphase_demo(option,line_style)
+%
+% option - 'fast'
+% 'real'
+%
+% line_style - just as plot, leave empty if you do not want lines between subbands
+%
+% Example:
+%
+% polyphase_demo('fast','r--');
+%
+% See also: reordering, partition, prefilter, polyphase_coefficients, polyphase
+%
+% (C) 2002 Astron, by M.van Veelen
+
+if nargin < 1
+ option='fast';
+end;
+
+figure(1) ;
+if strcmp(option,'fast') ; K1=16; L1=8 ; K2=16; L2=8 ; nms=50000; end ;
+if strcmp(option,'real') ; K1=64; L1=32 ; K2=64; L2=32 ; nms=1000000; end ;
+
+% x=chirp(t,50,1,150,'q');
+ dt=10E-6;
+ t=0:dt:dt*nms;
+ x=chirp(t,10,nms*dt,0.5/dt);
+% x=sin(2*pi*t*25) + sin(2*pi*t*50);
+
+% figure(1)
+% K=16;L=8;r=polyphase(x,K,L,K) ; subplot(121) ; imagesc(abs(r)') ; size(r)
+% K=16;L=8;r=polyphase(x,K,L,K,'shift') ; subplot(122) ; imagesc(abs(r)') ; size(r)
+% subplot(221) ; plot( t,x) ;
+
+
+y=polyphase(x,K1,L1,K1) ; %columns are subbands
+subplot(221) ; imagesc( abs(y)' ) ; ylabel( 'Subband after one stage of filtering') ;
+
+subband_demo=floor(K1*0.3);
+
+% plot one subband signal
+subplot(223) ; plot( real(y(:,subband_demo))) ; ylabel( 'Real part of a signal in one subband') ;
+axis( [ 0 length(y) min(min(x)) max(max(x))] )
+
+size(y)
+n_z = floor(length(y)/K2)-L2+1;
+z=zeros( n_z, K2*K1 );
+for k=1:K1
+ z_range=K2*(k-1)+[1:K2];
+ z(:,z_range)=polyphase(y(:,k),K2,L2,K2,'centre');
+end;
+
+subplot(222) ; imagesc( abs(z)' ) ; ylabel( 'Subband after the second stage of filtering') ;
+% plot sub band regions
+if nargin > 1
+ hold ;
+ for n=0:K1
+ plot( [ 1 length(z) ] , [ 1+n*K2 1+n*K2 ],line_style );
+ end;
+ hold ;
+end ;
+
+subplot(224) ; imagesc( abs(z(:,K2*(subband_demo)+[1:K2]) )') ; ylabel( 'The channels within this one subband') ;
+
+
diff --git a/SDP/SPP/MATLAB/polyphase_filter/polyphase_shifted.m b/SDP/SPP/MATLAB/polyphase_filter/polyphase_shifted.m
new file mode 100755
index 0000000000000000000000000000000000000000..8c41d3f55edc52e3d489ea091997dfcc61d3bdbf
--- /dev/null
+++ b/SDP/SPP/MATLAB/polyphase_filter/polyphase_shifted.m
@@ -0,0 +1,133 @@
+function [x,y,z]=polyphase_demo(option,line_style)
+%
+% polyphase_demo(option,line_style)
+%
+% option - 'fast'
+% 'real'
+%
+% line_style - just as plot, leave empty if you do not want lines between subbands
+%
+% Example:
+%
+% polyphase_demo('fast','r--');
+%
+% See also: reordering, partition, prefilter, polyphase_coefficients, polyphase
+%
+% (C) 2002 Astron, by M.van Veelen
+
+accuracy=2^(-14) ; dBmin=20*log10(accuracy) ;
+if nargin < 5 ; line_style = '' ; end ;
+if nargin < 1
+ option='fast';
+end;
+
+% figure(1) ;
+subplot(111) ; clf ;
+if strcmp(option,'fast') | strcmp(option,'fast-both');
+ K1=16; L1=8 ; M1=K1 ; N2=16 ; K2=N2; L2=8 ; M2=K2 ; nms=50000; dt=10^(-6);
+ x=zeros(1,nms); t=0:dt:dt*nms; x=chirp(t,10,nms*dt,0.5/dt);
+ [y,z]=polyphase_cascaded(x,K1,L1,M1,K2,L2,M2,N2);
+ size(y)
+ polyphase_demo_plots(x,y,z,dBmin,1,line_style)
+end ;
+
+if strcmp(option,'fast-nc') | strcmp(option,'fast-both');
+ K1=16; M1=8 ; L1=8 ; N2=16; K2=N2 * floor(K1/M1) ; L2=8 ; M2=K2 ; nms=50000; dt=10^(-6);
+ x=zeros(1,nms); t=0:dt:dt*nms; x=chirp(t,10,nms*dt,0.5/dt);
+ [y,z]=polyphase_cascaded(x,K1,L1,M1,K2,L2,M2,N2,'keep');
+ polyphase_demo_plots(x,y,z,dBmin,1,line_style)
+ [y,z]=polyphase_cascaded(x,K1,L1,M1,K2,L2,M2,N2,'cut');
+ polyphase_demo_plots(x,y,z,dBmin,2,line_style)
+end ;
+
+if strcmp(option,'real-nc') ;
+ K1=64; L1=32 ; M1=48; L1=floor(L1 * K1/M1) ; N2=64; K2=floor(N2*K1/M1/2)*2 ; L2=32 ;M2=K2 ; nms=1000000; dt=15*10^(-9);
+ x=zeros(1,nms); t=0:dt:dt*nms; x=chirp(t,10,nms*dt,0.5/dt);
+end;
+if strcmp(option,'real') ;
+ K1=32; L1=16 ; M1=K1 ; N2=64; K2=N2; L2=32 ; M2=K2 ; nms=640000; dt=15*10^(-9);
+ x=zeros(1,nms); t=0:dt:dt*nms; x=chirp(t,10,nms*dt,0.5/dt);
+end ;
+if strcmp(option,'hard') ;
+ K1=128; L1=32 ; M1=K1 ; N2=256; K2=N2; L2=32 ; M2=K2 ; nms=640000; dt=15*10^(-9);
+ x=zeros(1,nms); t=0:dt:dt*nms; x=chirp(t,10,nms*dt,0.5/dt);
+end ;
+
+% OLD x=chirp(t,50,1,150,'q');
+% OLD x=sin(2*pi*t*25) + sin(2*pi*t*50);
+
+
+% f1=256*64 ; f2=256*16;
+% x= [ sin(2*pi*t(1:length(t)/2-1)*f1) , sin(2*pi*t(length(t)/2:length(t))*f2) ] ;
+
+% figure(1)
+% K=16;L=8;r=polyphase(x,K,L,K) ; subplot(121) ; imagesc(abs(r)') ; size(r)
+% K=16;L=8;r=polyphase(x,K,L,K,'shift') ; subplot(122) ; imagesc(abs(r)') ; size(r)
+% subplot(221) ; plot( t,x) ;
+
+
+function [y,z]=polyphase_cascaded(x,K1,L1,M1,K2,L2,M2,N2,option)
+
+ if nargin<9 ; option='cut' ; end ;
+
+ y=polyphase(x,K1,L1,M1) ; y=y ./ max(max(abs(y))) ; %columns are subbands
+
+ if strcmp(option,'cut') ; dN2=floor((K2-N2)/2) ; else N2=K2; dN2=0; end;
+
+ n_z = floor(length(y)/K2)-L2+1;
+ z=zeros( n_z, N2*K1 );
+ for k=1:K1
+ z_range=N2*(k-1)+[1:N2];
+ size(1+dN2:K2-dN2) ;
+ size(z_range) ;
+
+ %r=polyphase(y(:,k),K2,L2,M2,'centre');
+ r=polyphase(y(:,k),K2,L2,M2);
+ %z(:,z_range)=r(:,1+dN2:K2-dN2);
+ z(:,z_range)=r(:,K2-dN2:-1:1+dN2);
+ end;
+ z=z./max(max(abs(z))) ;
+return ;
+
+function polyphase_demo_plots(x,y,z,dBmin,row_nr,line_style)
+ nr_columns = 3 ; nr_rows = 3 ;
+ if nargin < 5 ; row_nr = 1 ; end ;
+ K1=length(y')
+ K2=length(z')
+ subband_demo=1 ;% floor(K1*0.3);
+ % plot one subband signal
+ plot_nr=1+nr_columns*(row_nr-1) ; subplot(nr_rows,nr_columns,plot_nr) ;
+ plot( abs(y(:,subband_demo))) ; ylabel( 'Real part of a signal in one subband') ; view(0,90) ;
+ % axis( [ 0 length(y) min(min(x)) max(max(x))] )
+
+ plot_nr=2+nr_columns*(row_nr-1) ; subplot(nr_rows,nr_columns,plot_nr) ;
+ imagesc( max(dBmin, 20*log10(abs(y)')) ) ; ylabel( 'Subband after one stage of filtering') ; view(0,90) ;
+
+ plot_nr=3+nr_columns*(row_nr-1) ; subplot(nr_rows,nr_columns,plot_nr) ;
+ imagesc( max(dBmin, 20*log10(abs(z)') ) ) ;
+ ylabel( 'Subband after the second stage of filtering') ; view(0,90) ;
+ % axis( [ 0 length(y) min(min(z)) max(max(z))] )
+ % plot sub band regions
+ if ~strcmp(line_style,'') ; hold ;
+ for n=0:K1
+ plot( [ 1 length(z) ] , [ 1+n*K2 1+n*K2 ],line_style );
+ end; hold ;
+ end ;
+
+% subplot(344) ; mesh( max(dBmin,20*log10(abs(z(:,K2*(subband_demo)+[1:K2])))')) ; view(0,90) ;
+% ylabel( 'The channels within this one subband') ; %
+
+ colormap gray;
+
+return ;
+
+if 0
+ X1=fft(x(partition(length(x),K1)')) ; X1=X1./max(max(abs(X1)));
+ subplot(333); mesh(max(dBmin,20*log10(abs(X1)))) ; view(0,90) ;
+ ylabel( 'FFT result') ;
+
+ X2=fft(x(partition(length(x),K1*K2)')) ; X2=X2./max(max(abs(X2)));
+ subplot(339); mesh(max(dBmin,20*log10(abs(X2)))) ; view(0,90) ;
+ ylabel( 'FFT result') ;
+end;
+
\ No newline at end of file
diff --git a/SDP/SPP/MATLAB/polyphase_filter/prefilter.m b/SDP/SPP/MATLAB/polyphase_filter/prefilter.m
new file mode 100755
index 0000000000000000000000000000000000000000..5c191405910ddc3cd0ab09e0ba6896a7ec31b25a
--- /dev/null
+++ b/SDP/SPP/MATLAB/polyphase_filter/prefilter.m
@@ -0,0 +1,53 @@
+function [ym] = prefilter(x,K,L,M,centre,h)
+%
+% function [Y] = prefilter(x,K,L,M,centre,h)
+%
+% K - the number of channels
+% L - filter length
+% M - the downsampling rate
+% centre - if centre='yes' the output is shifted to the centre (default centre='no')
+% h - replaces the default fircls1 filter window with the one specified
+%
+% Compute the polyphase filter coefficients p and applies them the the data x, i.e. the
+% data is decimated by a factor M and divided into K channels. The prototype filter can
+% be supplied through h, by default fircls1 is used to derive the prototype filter window.
+% The coefficients care derived from this prototype windows and convolved with x, hence
+% a number of output vectors are computed which can be transformed into the real filter
+% outputs using an DFT. The polyphase filter structure comprises of two stages:
+%
+% STAGE 1: Convolution
+% STAGE 2: Fourier Transform
+%
+% The first stage is implemented in the function. The output is an N X K matrix with a
+% channel in each columm, providing K columns; the length N, is floor((length(x))/K)-L+1,
+% in case no case the data is critically sampled (K=M).
+%
+% See also: reordering, partition, prefilter, polyphase, polyphase_coefficients
+%
+% (C) 2002 Astron, by M.van Veelen
+
+% Settings for finite word-length coding
+
+% INITIALIZE: compute the filter prototype coefficients of the
+
+if nargin < 5 ; centre='no' ; end ;
+
+if nargin == 6
+ [i,h]=polyphase_coefficients(K,L,M,'indices',centre,h);
+else
+ [i,h]=polyphase_coefficients(K,L,M,'indices',centre);
+end;
+
+n_y = floor((length(x)) / K )-L+1;
+parts=partition(length(x),K,0)';
+xd=x(parts) ;
+y = zeros(1, K);
+ym= zeros(K, n_y);
+
+% STAGE 1: apply the pre-filter
+for m = 1 : n_y
+ y=1/L * sum(( xd(i+(m-1)*K) .* h(i) )' );
+ % When p=h(i) this is equivalent to y=1/L * sum(( xd(:, m : m + L - 1) .* p )' );
+ % y = double(uencode(y, quant_inputfft, 1, 'signed')) / 2^(quant_inputfft - 1); % Quantise the input of the FFT
+ ym(:,m)=y';
+end ;
\ No newline at end of file
diff --git a/SDP/SPP/MATLAB/polyphase_filter/reordering.m b/SDP/SPP/MATLAB/polyphase_filter/reordering.m
new file mode 100755
index 0000000000000000000000000000000000000000..4bbed0027a7a185109b161daab15ac2b4ddf1a63
--- /dev/null
+++ b/SDP/SPP/MATLAB/polyphase_filter/reordering.m
@@ -0,0 +1,79 @@
+function [y] = reordering( x,M,L,option )
+%
+% function [Y] = reordering( X,M,L,OPTION )
+%
+% input:
+% X - input vector
+% M - number of channels and downsample rate
+% L - upsample rate [default 1]
+% OPTION - 'block'
+% 'streaming'
+% 'demo'
+% 'test'
+%
+% output:
+% Y - reordered data matrix in case X is a vector:
+%
+% Y has M rows and L*length(x)/M columns
+%
+% The input signal is decimated by a factor M, thus obtaining M different channels
+% with a sampling frequency of fs/M when fs is the sampling frequency of the input.
+% In the default case L=1 the samples are distributed uniformly, for higher values
+% of L each channel is upsampled by L, thus within each channel the samples are
+% M/L apart, e.g. with reordering([1:30],8,1.5) the returned matrix will be:
+%
+% 1 6 11 16 21
+% 2 7 12 17 22
+% 3 8 13 18 23
+% 4 9 14 19 24
+% 5 10 15 20 25
+% 6 11 16 21 26
+% 7 12 17 22 27
+% 8 13 18 23 28
+%
+% This funtion is implemented using partition. The indices of the input matrix at the
+% location of the output matrix of the reordering function can be obtained with:
+% partition(length(x),M,floor(M*(L-1))) , this indices can also be used to obtain Y from X:
+%
+% reordering(x,M,L) == x(partition(length(x),M,floor(M*(L-1))))'
+%
+% See also: partition
+%
+% (C) 2002 Astron, by M.van Veelen
+
+if nargin < 4 ;
+ option = 'block' ;
+end ;
+
+if nargin < 3
+ L=1;
+end;
+
+if strcmp(option,'streaming')
+ error('option streaming is not implemented yet');
+end ;
+
+if strcmp(option,'demo')
+ display( 'x=1:100;y=reordering(x,1.3,1);mesh(y)' ) ;
+ x=1:100;y=reordering(x,16,1.25);mesh(y);
+ return ;
+end ;
+
+if strcmp(option,'test')
+ if reordering(x,M,L) == x(partition(length(x),M,floor(M*(L-1))))'
+ display('TEST SUCCESSFUL : ASSERT{reordering(x,M,L)== x(partition(length(x),M,floor(M*(L-1))))''}') ;
+ end ;
+ return ;
+end ;
+
+n_x = length(x) ;
+indices=partition( n_x, M, floor(M*(L-1)),'range' );
+n_y=length(indices');
+y=zeros(M,n_y) ;
+
+for i=1:n_y
+ y(1:M,i) = x(indices(i,1):indices(i,2))';
+% DEBUG t=x(indices(i,1):indices(i,2))
+% DEBUG y(1:M,i) = t';
+end ;
+