Cross-Spectral Analysis of Two Signals in MATLAB

% comments on the code:/n/n% The code is used to analyze two signals (data1 and data2) using cross-spectral/n% analysis. The signals are first filtered using a Butterworth filter/n% within a specified frequency range. The filtered signals are then/n% divided into windows of length ts and a Hanning window is applied to each/n% window. The Fourier transform is then applied to each window and the/n% cross-spectral analysis is performed./n/n% The output of the code includes the power spectra of both signals, the/n% cross-power spectrum, the coherence spectrum, and the phase spectrum./n% These spectra are plotted against frequency./n/n% The code uses built-in functions from MATLAB for filtering, Fourier/n% transform, and windowing. It also uses detrend function to reduce the/n% ripple of data./n/n% The code can be used for various applications, such as analyzing EEG or/n% ECG signals, seismic signals, or any other signals that can be/n% represented in the time domain./ns = 4000; % the width of filter windows, ts numbers per windows/nts = 4000; % the width of filter windows, ts numbers per windows/nwx = hanning(ts); % generates a hanning filter coefficients/nfs = 2000; % the sampling rate of probes, [kHz]/nfre_cal = fix((length(data1)-ts)/(ts/2)); % calculate the frequency of analysis/nrange_butter = [1, 999999]/1000000; % the frequency range of butterworth filter/n%% butterworth filter/n[coef1, coef2] = butter(1, range_butter); % generates the coefficients of butterworth filter/ndata1_but = filter(coef1, coef2, filter(coef1, coef2, data1)); % data processed by butterworth filter/ndata2_but = filter(coef1, coef2, filter(coef1, coef2, data2)); % data processed by butterworth filter/n%% FFT/ndata1_fft = zeros(fre_cal, ts); % generates a box for the data after FFT/ndata2_fft = zeros(fre_cal, ts);/ndisp('系综平均数：')/ndisp(fre_cal)/nfor i = 1 : fre_cal/n data1_fft(i, :) = data1_but((i-1)*ts/2 + 1 : (i-1)*ts/2 + ts); % read the data for every filter windows/n data2_fft(i, :) = data2_but((i-1)*ts/2 + 1 : (i-1)*ts/2 + ts);/n data1_detr(i, :) = detrend(data1_fft(i, :)); % reduce the ripple of data/n data2_detr(i, :) = detrend(data2_fft(i, :));/n data1_fft(i, :) = fft(data1_detr(i, :).wx', ts);/n data2_fft(i, :) = fft(data2_detr(i, :).wx', ts);/nend/n% calculate the frequency series. % f = nfs/N/n% N is the number of data points, n is a vector (0 to N - 1)/n% fs is sampling rates/nfrequency = (0 : (ts/2 - 1))(fs/(ts - 2));/n%% calculate/nmx12 = sum(abs(data1_fft(:, 1:ts/2)).*abs(data2_fft(:, 1:ts/2)))/fre_cal; % calculate (Σ|x1|·|x2|)/fre_cal/n% figure, plot(frequency(2 : ts/2), mx12(2:end)), xlim([fBegin, fEnd])/nphase = angle(sum(data1_fft(:, 1:ts/2).*conj(data2_fft(:, 1:ts/2))))/pi; % calculate the phase of all filter windows, [rad]/npower_cross = abs(sum(data1_fft(:, 1:ts/2).*conj(data2_fft(:, 1:ts/2))))/fre_cal; % calculate the average cross power of all filter windows/npower1 = sum(((abs(data1_fft(:, 1:ts/2))).^2))/fre_cal; % calculate the power of data1/npower2 = sum(((abs(data2_fft(:, 1:ts/2))).^2))/fre_cal; % calculate the power of data2/ncorherency = power_cross./mx12; % calculate corherency/n% fre_cal/n% length(phase)/nfreOutput = frequency(2:ts/2);/nphaseOutput = phase(2 : end);/npower1Output = power1(2:end);/npower2Output = power2(2:end);/npowerCsOutput = power_cross(2:end);/nCorrOutput = corherency(2:end);/n%% figures/nfigure/nsubplot(511)/nplot(frequency(2 : ts/2), smooth(power1(2 : end), 3), 'LineWidth', 1), xlim([fBegin, fEnd])/nylabel('power1 (a.u)')/nsubplot(512)/nplot(frequency(2 : ts/2), smooth(power2(2 : end), 3), 'LineWidth', 1), xlim([fBegin, fEnd])/nylabel('power2 (a.u)')/nsubplot(513)/nplot(frequency(2 : ts/2), power_cross(2 : end), 'LineWidth', 1), xlim([fBegin, fEnd])/nylabel('cross power (a.u)')/nsubplot(514)/nplot(frequency(2 : ts/2), smooth(corherency(2 : end), 5), 'LineWidth', 1), xlim([fBegin, fEnd])/nylabel('corherency')/nsubplot(515)/nplot(frequency(2 : ts/2), smooth(phase(2 : end), 3), 'LineWidth', 1), xlim([fBegin, fEnd])/nhold on/nplot([0, 500], [0, 0], 'b--', 'linewidth', 1), ylim([-1, 1])/nylabel('phase (rad//pi)')/nxlabel('f (kHz)')