Skip to content
Snippets Groups Projects
Commit 3c54a03d authored by gerdes's avatar gerdes
Browse files

This commit was generated by cvs2svn to compensate for changes in r693, 

which included commits to RCS files with non-trunk default branches.
parent 776ee075
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 693 additions and 0 deletions
SDP/SPP/MATLAB/polyphase_filter/partition.m 0 → 100755
+ 54
− 0
View file @ 3c54a03d
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;
SDP/SPP/MATLAB/polyphase_filter/polyphase.m 0 → 100755
+ 59
− 0
View file @ 3c54a03d
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);
SDP/SPP/MATLAB/polyphase_filter/polyphase_coefficients.m 0 → 100755
+ 112
− 0
View file @ 3c54a03d
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);
SDP/SPP/MATLAB/polyphase_filter/polyphase_compare.m 0 → 100755
+ 136
− 0
View file @ 3c54a03d
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 ;
SDP/SPP/MATLAB/polyphase_filter/polyphase_demo.m 0 → 100755
+ 67
− 0
View file @ 3c54a03d
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') ;
SDP/SPP/MATLAB/polyphase_filter/polyphase_shifted.m 0 → 100755
+ 133
− 0
View file @ 3c54a03d
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
SDP/SPP/MATLAB/polyphase_filter/prefilter.m 0 → 100755
+ 53
− 0
View file @ 3c54a03d
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
SDP/SPP/MATLAB/polyphase_filter/reordering.m 0 → 100755
+ 79
− 0
View file @ 3c54a03d
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 ;
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment