Muscle synergy analysis is a kind of modularized decomposition of muscles during exercise controlled by the central nervous system (CNS). It can not only extract the synergistic muscles in exercise, but also obtain the activation states of muscles to reflect the coordination and control relationship between muscles. However, previous studies have mainly focused on the time-domain synergy without considering the frequency-specific characteristics within synergy structures. Therefore, this study proposes a novel method, named time-frequency non-negative matrix factorization (TF-NMF), to explore the time-varying regularity of muscle synergy characteristics of multi-channel surface electromyogram (sEMG) signals at different frequency bands. In this method, the wavelet packet transform (WPT) is used to transform the time-scale signals into time-frequency dimension. Then, the NMF method is calculated in each time-frequency window to extract the synergy modules. Finally, this method is used to analyze the sEMG signals recorded from 8 muscles during the conversion between wrist flexion (WF stage) and wrist extension (WE stage) movements in 12 healthy people. The experimental results show that the number of synergy modules in wrist flexion transmission to wrist extension (Motion Conversion, MC stage) is more than that in the WF stage and WE stage. Furthermore, the number of flexor and extensor muscle synergies in the frequency band of 0–125 Hz during the MC stage is more than that in the frequency band of 125–250 Hz. Further analysis shows that the flexion muscle synergies mostly exist in the frequency band of 140.625–156.25 Hz during the WF stage, and the extension muscle synergies appear in the frequency band of 125–156.25 Hz during the WE stage. These results can help to better understand the time-frequency features of muscle synergy, and expand study perspective related to motor control in nervous system.