We propose a novel approach for demodulating continuous phase modulation (CPM) signals based on the modeling of the instantaneous phase as a piecewise polynomial-phase function. The polynomial modeling can be a good approximation for currently used modulations or it can be exact if the shaping pulse is chosen to be a piecewise polynomial function. The crucial step in the demodulation process is then the estimation of the polynomial coefficients, which is carried out using the so called product high order ambiguity function (PHAF). The proposed approach is suboptimal with respect to the optimal maximum likelihood sequence estimation (MLSE) method, but is much simpler to implement and offers important advantages such as independence of initial phase, tolerance to Doppler shift, and time-offset, blind channel identification. We show theoretical results concerning the minimum distance among sequences, which leads to a lower bound on the error probability, together with some simulation results.