A fully coupled thermoelastic framework is formulated to cope with the free vibration response of anisotropic multilayered plates in three dimensions. The laminated structure consists of homogeneous laminae of arbitrary thickness and width under simply supported edge conditions in thermal environment. The general and exact field expressions of the temperature, heat flux, displacement and stress components are expressed in terms of double Fourier series expansions in any rectangular plate, which lead to the extended Stroh formalism with thermomechanical coupling effects in a concise and compact matrix form. Different imperfect interface conditions are introduced to characterize specific structural and thermal contact properties at the bounding interfaces, and further to determine the finite complex valued coefficients in the suitable series relations. The complete time-harmonic solutions in the laminated composites in the presence of perfect/imperfect interfaces are recursively obtained by means of the modified dual variable and position technique with explicit layer-to-layer transfer matrices. Results are obtained for different layups, length-to-thickness ratios and interfacial boundary conditions for two application examples, namely the graphite/epoxy cross-ply composites and the thermal barrier coatings on superalloys, without suffering from numerical exponential instability. These investigations reveal that the natural frequencies and first and higher vibration mode shapes of the multilayered structures can be considerably affected by increasing the environmental temperature and the severity of the interfacial imperfections. Since the through-thickness stress distribution in 2, 5, and 10 layered composites appears to be strongly correlated to the layups, such modal stress analysis could be exploited to locate the fatigue hotspots operated in dynamic structures and to guide the structural design of aircraft and spacecraft composite laminates subjected to residual vibrations.