A method for calculating energy levels and wavefunctions of bound states associated with deep impurity centres in semiconductors is presented. The impurity system is replaced by a crystal, whose unit cell consists of the impurity surrounded by a large number of host crystal atoms. The impurity bands for this crystal, represented by a real potential, are calculated by expanding in orthogonalised Bloch waves. To find the discrete level and the corresponding isolated impurity wavefunctions, a tight-binding extrapolation scheme has been developed. As a first application, the method is applied to substitutional neutral zinc in silicon, using 54 atoms per unit cell. The total impurity bandwidth is found to be quite large, 1.5 eV, arising from a discrete level at about 0.2 eV above the top of the valence band. The many-band nature and probability distribution of the impurity of Bloch states are explored, showing moderate localisation in both real and reciprocal space.