A general formulation of the method of the reverberation-ray matrix (MRRM) based on the state space formalism and plane wave expansion technique is presented for the analysis of guided waves in multilayered piezoelectric structures. Each layer of the structure is made of an arbitrarily anisotropic piezoelectric material. Since the state equation of each layer is derived from the three-dimensional theory of linear piezoelectricity, all wave modes are included in the formulation. Within the framework of the MRRM, the phase relation is properly established by excluding exponentially growing functions, while the scattering relation is also appropriately set up by avoiding matrix inversion operation. Consequently, the present MRRM is unconditionally numerically stable and free from computational limitations to the total number of layers, the thickness of individual layers, and the frequency range. Numerical examples are given to illustrate the good performance of the proposed formulation for the analysis of the dispersion characteristic of waves in layered piezoelectric structures. multilayered piezoelectric structures, the method of reverberation-ray matrix, guided waves, dispersion characteristic Piezoelectric materials have been applied in many important fields such as geophysics, electronics, communication, instrumentation and nondestructive evaluation and testing of materials [1] . Numerous ultrasonic devices, based on either surface acoustic wave (SAW) or bulk acoustic wave (BAW), have been developed using piezoelectric materials for a variety of applications [2,3] . These devices are usually fabricated with a layered geometry, because of its high sensitivity, great bandwidth, and enhanced reception characteristics [4][5][6][7][8][9][10][11] . Currently, in order to satisfy the increasing demand of higher frequency, higher performance, smaller size, lower cost and lower energy consumption, thin film piezoelectric stack structures are widely used with the maturation of thin film technologies. Thin film bulk acoustic resonator (TFBAR) and solidly mounted resonator (SMR) are two typical examples [12][13][14] . For the sake of design and optimization, appropriate theoretical models and efficient solution methods are desired to investigate the behaviors of wave propagation in these devices. It is common to use layered model consisting of piezoelectric and non-piezoelectric layers stacked in a certain sequence for these piezoelectric devices [4][5][6][7][8][9][10][11][12][13] . So far, various matrix analysis methods have been proposed in the study of wave propagation in multilayered piezoelectric structures [1][2][3][4][5][6][7][8][9][10][11][12][15][16][17][18][19][20][21][22][23][24][25][26][27] , such as the finite element method (FEM) [4,10] , the transfer matrix method (TMM) [1][2][3]15,16] , the scattering-matrix method [8] , the surface impedance