fftd(t)

function [b, d] = fftd(t)
%Dimensionless Discrete Fourier Transform
%
%Usage: [B, D] = fftd(T)
%
%    FFTD(T) performs a "Dimensionless DFT" on the columns of T.
%    T may be real or complex. T may be any size.
%    Returns the dimensionless (unitless) vector D
%    and the universal Basis matrix B such that
%
%    Time = B*D,   and  Freq = Inv(B)*D
%
%    Where Freq = fft(Time)/n is the conventional DFT pair.
%    Note: B has already been scaled such that mean(t) = f(1)
%
%    The Basis matrix, B, is always the same for the same sized T
%    with Time = B*B*Freq
%
%    D is the dimensionless primitive representation of
%    power signal T and includes Time and Frequency
%    information simultaneously.
%
%    If T is odd the Time=kt*Real(D) and Freq=kf*Imag(D)
%
%    Parseval's Theorem still holds (with correct scaling) for D
%
%    Tested under version 6.0.0.88 (R12)
%    See also FFT, IFFT, FFTH, FFT2, IFFT2, FTSHIFT

% Paul Godfrey
% [e-mail]
% 6-6-2003

%revision history
%6-6-2003 fixed bug in computing b
%6-6-2003 allowed for complex t of any size

sizeoft=size(t);
n=sizeoft(1);
t=t(:);
t=reshape(t,n,length(t)/n);

m=fft(eye(n))/n;
bi=sqrtm(m);
b=inv(bi);
d=bi*t;
d=reshape(d(:),sizeoft);

%we have  t =  b *columns_of_d
%and      f =  bi*columns_of_d

return

%B is just the matrix square root of the unit circle
plot(b*b,'*');axis image;grid on

% A demo of this function is:
clc
clear all
close all
warning off
format short g
big=2^45;

n=64; 
k=(0:n-1).';
w=k*2*pi/n;

t=0;
for in=1:2:11;
    t=t+sin(in*w)/in;
end
%for in=13:2:23;
%    t=t+cos(in*w)/in;
%end

[b,d]=fftd(t);
f=b'*d/sqrt(n);

f=round(f*big)/big;
d=round(d*big)/big;
t=round(t*big)/big;
%sin is an odd function so everything separates
fpt=[8*f 4*d t];

figure(1);plot3(k,imag(t),abs(t).*sign(t),'b');hold on;
figure(1);plot3(k,imag(d),abs(d).*sign(d)*4,'r*');axis tight;
grid on;rotate3d;view([0,0]);

figure(2);plot3(k,imag(f),abs(f).*sign(f),'b');hold on;
figure(2);plot3(k,imag(d)/2,abs(d).*sign(d),'r*');axis tight;
grid on;rotate3d;view([0,90]);

figure(3);plot3(k,imag(d),abs(d).*sign(d),'r*');axis tight;
grid on;rotate3d;view([-15,45]);
disp('Use the mouse to examine/move the overlaid data');

clear all
n=8;
t=rand(n,1);
f=fft(t)/n;
[b,d]=fftd(t);
tt=[t b*d]
ff=[f b'*d/sqrt(n)]

figure(3)
set(gcf,'doublebuffer','on')
for el=0:5:90
    view(0,el);
    drawnow
end
for el=90:-5:0
    view(0,el);
    drawnow
end

return