In this study, our goal is to investigate the multivariate rhythmic characteristics of electroencephalogram (EEG) and electrooculogram (EOG) data in adaptive frequency scales for drowsiness detection. The empirical Fourier decomposition (EFD) method has been extended to multivariate signals to eliminate the mode alignment issue in multichannel signals. The proposed multivariate EFD (MEFD) is designed to adapt mode alignment to extract mutual characteristics across multi-channel and delivered channel-aligned modes. This novel technique utilizing an adaptive zero-phase filter bank, optimized boundary estimation, and a MEFD-based multi-model framework has been introduced to enhance drowsiness recognition performance by addressing invariant sampling rates across EEG and EOG modalities. In order to test effectiveness, the proposed MEFD has been investigated using multichannel synthetic signals, and EEG and EOG signals from the SEED-VI drowsiness database. Further feature fusion matrix (FFM) has been formed to extract mutual information from multivariate multimodal EEG and EOG signals. Finally, the obtained FFM features have been tested using six well-known classifiers. With 5-fold cross-validation, our proposed MEFDbased classification framework has demonstrated exceptional performance in drowsiness state classification, achieving the highest accuracy of 82.3% when compared to state-of-the-art methods.