This work presents a fast Fourier transform (FFT) based method that can be used to model interface decohesion. The inability of an FFT solver to deal with sharp interfaces discards the use of conventional cohesive zones to model the interfacial mechanical behaviour within this framework. This limitation is overcome by approximating sharp interfaces (e.g. grain/phase boundaries) with an interphase. Within the interphase, the background plastic constitutive behaviour is inherited from the respective adjacent grains. The anisotropic kinematics of the decohesion process is modelled using a damage deformation gradient that is constructed by mapping the opening strains (in normal and tangential modes) to the associated projection tensors. The degradation (damage) of the interfacial opening resistances is modelled via a dimensionless nonlocal damage variable that prevents localisation of damage within the interphase. This nonlocal variable results from the solution of a gradient damage based regularisation equation within the interphase subdomain. The damage field is constrained to the interphase region by applying a relatively large penalisation on the damage gradients inside the interphase. The extent of nonlocality ensures that the damage is largely uniform in the direction perpendicular to the interphase, thus making its thickness the theoretical lengthscale for dissipation. To achieve model predictions that are objective with respect to the interphase thickness, scaling relations of the model parameters are proposed. The numerical performance is shown for a uniform opening case and then for a propagating crack. Finally, an application to an artificial polycrystal is shown.