How to set up the beginning of the coordinate system in the bottom left corner for Wigner-Hough Transform?

Hello! I'm trying to calculate Wigner-Hough Transform. I saw some examples, especially the "htl.m" file in the Time-Frequency Toolbox. But the beginning of the coordinate system is set up in the middle of the image and I don't understand why. I'd like to set it in the bottom left corner. Obviously, I'm missing something. I wrote this with the help of "htl.m", but it doesn't give me the correct result. The function wvch is also provided after the main program. Any suggestions will be appreciated!

clear t1=0:0.001:1; N1 = length(t1); sig = chirp(t1,0,0.8,150); sig = hilbert(sig); % Signal in time domain figure(1); subplot(221); plot(t1,real(sig)); xlabel('t'); ylabel('A'); % Wigner-Ville Distribution [tfr,t,f] = wvch(sig); f = f * (N1-1)/2/N1; t = t * 1/1000; [F, T] = meshgrid(f, t); subplot(223); mesh(F, T, abs(tfr)); subplot(224); contour(f,t,abs(tfr)); xlabel('f'); ylabel('t');

[Xmax,Ymax] = size(tfr); rhomax=sqrt(Xmax^2+Ymax^2)/2; M = Xmax; N = Ymax; deltar=rhomax/(M-1); deltat=2*pi/N;

ht=zeros(M,N); Max=max(max(tfr)); for x = 0:Xmax-1, for y = 0:Ymax-1, if abs(tfr(x+1,y+1))>Max/20, for theta = 0:deltat:(2*pi-deltat), rho = x*cos(theta) - y*sin(theta) if ((rho>=0)&&(rho<=rhomax)), ht(round(rho/deltar)+1,round(theta/deltat)+1) =... ht(round(rho/deltar)+1,round(theta/deltat)+1)+... tfr(x+1,y+1); end end end end end rho=0:deltar:rhomax; theta=0:deltat:(2*pi-deltat);

figure; mesh(rho, theta*180/pi, abs(ht)); xlabel('rho'); ylabel('theta'); figure; image(rho,theta*(180/pi), abs(ht)); xlabel('rho'); ylabel('theta');

%%%%%% wvch %%%%%% function [tfr, t, f] = wvch(x, t, N) %WVCH Wigner-Ville distribution from a Chinese book % [TFR,T,F]=WVCH(X,T,N) computes the Wigner-Ville distribution % of a discrete-time signal X, which is given in the form (1xN), % or the cross Wigner-Ville representation between two signals. % % X : signal if auto-WV, or [X1,X2] if cross-WV. % T : time instant(s) (default : 1:length(X)). % N : number of frequency bins (default : length(X)). % TFR : time-frequency representation. % F : vector of normalized frequencies.

if (nargin == 0), error('At least one parameter required'); end;

[xrow,xcol] = size(x); if xrow ~= 1, error('X must have one row, (1xN) form'); end;

if (nargin == 1), t=1:xcol; N=xcol ; elseif (nargin == 2), N=xcol ; end;

if (N<0), error('N must be greater than zero'); end;

[trow,tcol] = size(t); if (trow~=1), error('T must only have one row'); end;

N1 = length(x)+rem(N,2); length_FFT = N1; % Take an even value of N R = zeros(N,N); for n = 0:N-1 M = min(n, N-1-n); for k = 0:M R(n+1, k+1) = x(n+k+1) * conj(x(n-k+1)); end for k = N-1 : -1 : N-M R(n+1, k+1) = conj(R(n+1, N-k+1)); end end tfr = zeros(N1, N1); for n = 0: N-1 temp = fft(R(n+1,:), length_FFT)/(2*length_FFT); temp = fftshift(temp); tfr(n+1, :) = temp; end t1 = linspace(0, 1, N1); f = (t1*(N-1) - N1/2); t = (0 : N1-1);