examp_pgauss.m

linear time－frequency toolbox

%EXAMP_PGAUSS  How to use PGAUSS
%
%  This script illustrates various properties of the
%  Gaussian function.
%
%  FIGURE 1 Window+Dual+Tight
%
%   This figure shows an optimally centered Gaussian for a 
%   given Gabor system, its canonical dual and tight windows
%   and the DFTs of these windows.
%
%  FIGURE 2 Framebounds sweep.
%
%   This figure shows how the framebounds behave for a Gabor frame
%   with a fixed lattice, but with Gaussian windows of different
%   widths.
%

disp('Type "help examp_pgauss" to see a description of how this example works.');

% A quick test: If the second input parameter to
% pgauss is not specified, the output will be
% invariant under an unitary DFT. Matlabs FFT is does not
% preserve the norm, so it must be scaled a bit.

L=128;
g=pgauss(L);

disp('');
disp('Test of DFT invariance: Should be close to zero.');
norm(g-dft(g))

% Setup parameters and length of signal.
% Note that it must hold that L=M*b=N*a for some integers
% b and N, and that a <= M
L=72;  % Length of signal.
a=6;   % Time shift.
M=9;   % Number of modulations.

% Calculate the frequency shift.
b=L/M;

% For this Gabor system, the optimally concentrated Gaussian
% is given by
g=pgauss(L,a/b);

% This is not invarient with respect to a DFT, but it is still
% real and whole point even
disp('');
disp('The function is WP even. The following should be 1.');
iseven(g)

disp('Therefore, its DFT is real.');
disp('The norm of the imaginary part should be close to zero.');
norm(imag(dft(g)))

% Calculate the canonical dual.
gdual=candual(g,a,M);

% Calculate the canonical tight window.
gtight=cantight(g,a,M);

% Plot them:

% Standard note on plotting:
%
% - All windows have real DFTs, but Matlab does not
%   always recoqnize this, so we have to filter away
%   the small imaginary part by calling REAL(...)
%
% - The windows are all centered around zero, but this
%   is not visually pleasing, so the window must be
%   shifted to the middle by an FFTSHIFT
%

gf_plot      = fftshift(real(dft(g)));
gdual_plot   = fftshift(gdual);
gdualf_plot  = fftshift(real(dft(gdual)));
gtight_plot  = fftshift(gtight);
gtightf_plot = fftshift(real(dft(gtight)));
figure(1);

subplot(3,2,1);
x=(1:L).';
plot(x,fftshift(g),'-',...
     x,circshift(fftshift(g),a),'-',...
     x,circshift(fftshift(g),-a),'-');
title('g=pgauss(72,6/8)');
legend('off');

subplot(3,2,2);
plot(gf_plot);
title('g, frequency domain');
legend('off');

subplot(3,2,3);
plot(gdual_plot);
title('Dual window of g');
legend('off');

subplot(3,2,4);
plot(gdualf_plot);
title('dual window, frequency domain');
legend('off');

subplot(3,2,5);
plot(gtight_plot);
title('Tight window generated from g');
legend('off');

subplot(3,2,6);
plot(gtightf_plot);
title('tight window, frequency domain');
legend('off');

%
% Plot the behaviour of the framebounds when the parameter
% of pgauss is swept over a range of values. 

% Number of plotting points.
npoints=100;

% The range of parameter
w=logspace(-1,1,npoints)*a/b;

% Calculate the framebounds for each plotting point.
ar=zeros(npoints,1);
br=zeros(npoints,1);
for ii=1:npoints
  [ar(ii),br(ii)]=gfbounds(pgauss(L,w(ii)),a,M);
end;


figure(2);

loglog(w,ar,'-',w,br,'-');

title('Framebounds');
xlabel('Parameter for PGAUSS');
legend('off');