We introduce a new method of deriving numerical algorithms for Linear Canonical Transforms (LCT) based on matrices, which act on phase space and distort the shape of the Wigner Distribution Function. Special cases of the LCT include the Fourier Transform (FT), the fractional Fourier Transform (FRT), the Fresnel Transform (FST). We show that many of the existing algorithms, which have been discussed in the literature, can be derived efficiently using this method. They can also be optimised and the relationship between them is discussed. In the case of the FRT all of the existing algorithms can be made index additive and reversible using correct amounts of interpolation and decimation. We derive many new algorithms for the LCT and show the means for deriving many more.