The real and imaginary parts of the Fourier transform of a causal signal can be estimated from one another using the Hilbert transform. The numerical computation can be carried out by a fast Fourier Transform or the convolution of the data with an appropriate Hilbert kernel. However, magnetotelluric data are unequally sampled because logarithmic frequency axes are conventionally used. A method was developed for the Hilbert transformation of unequally sampled data and applied to the real and imaginary parts of the frequency normalized impedance to examine if the dispersion relations are satisfied for a given impedance tensor. A version of the Shannon sampling theory was generalized for the unequally sampled data. The method involves the reconstruction of the measured data by a linear combination of an analytic interpolation function distributed along the horizontal axis. Then, the decomposition coefficients that provide a fit between the measured and reconstructed data are solved using the linear least-squares method with singular value decomposition. Finally, the discrete Hilbert transformed data are constructed at any desired abscissa values by the sum of the analytic Hilbert transform of the interpolation function using the decomposition coefficients solved in the previous step. Examples from the magnetotelluric synthetic and field data, including out-of-quadrant phase cases, are given. In principle, the method can also be applied to other integral transforms.