How to find unwrapped phase of a cross-spectral at frequencies of 0.1 0.2 0.3 ... 0.9 Hz?

3 views (last 30 days)
Hello everyone. I am new to matlab and also to the signal processing. I have to find the unwrapped phase of the cross-spectral. I have two signals with the sampling frequency of 100Hz. I have taken the fft of the signals and cross correlated them in frequency domain to make cross-spectral. Now the thing I am looking for is to find the unwrapped phase of the cross-spectral and plot its values against frequencies of 0.1 0.2 0.3 upto 0.9 Hz. The first colum of the text files of the signals is time and the second column is for the corresponding amplitude. Can someone guide me about that?
sig1 = readtable("total.txt");
sig2 = readtable("HAK_PHA_2017_11_26.txt");
s1 = table2array(sig1);
s2 = table2array(sig2);
L1 = length(s1);
L2 = length(s2);
F1 = fft(s1(:,2));
F2 = fft(s2(:,2));
Cross_spectral = corr(F1,F2)

Sign in to comment.

Accepted Answer

Mathieu NOE
Mathieu NOE on 6 Sep 2022
minor hint : readmatrix suffices to read numerical txt files ; readtable and conversion from table to array is unnecessary and overkill
second topic : you want a result with a frequency spacing of 0.1 Hz; as the sampling rate is 100 Hz, this means we need at least 1000 samples ; need zero padding as your file contains only 601 samples
tried your code vs my own method and both results match (fortunately !)
here we are :
s1 = readmatrix("total_1770_1776.txt");
s2 = readmatrix("HAK_PHA_2017_11_26_1770_1776.txt");
Samples = length(s1);
Ts = mean(diff(s1(:,1))); % Sampling Timeperiod
Fs = 1/Ts; % Sampling frequency
% padd with zeros to have at least 1000 samples so that we can achieve df =
% 0.1 Hz (df = Fs / nfftn and nfft is maximised , equal to number of
% samples available);
samples = 1000;
sig1 = zeros(samples,1);
sig1(1:length(s1)) = s1(:,2);
sig2 = zeros(samples,1);
sig2(1:length(s2)) = s2(:,2);
%% method 1 (yours)
F1 = fft(sig1);
F2 = fft(sig2);
Corrs1s2 = (ifft(F1.*conj(F2))); % Cross-Correlation
CS = fft(Corrs1s2); % Cross-spectrum
%% method 2 (mine)
% FFT parameters
nfft = length(sig1);
% window = hanning(nfft);
window = ones(nfft,1); % box window
noverlap = 0;
[Txy,Cxy,Pxy,f] = tfe_func(sig2,sig1,nfft,Fs,window,noverlap);
% as df = 0.1 Hz , select the first 10 values (frequencies of 0 0.1 0.2 0.3 ... 0.9 Hz)
ind = 1:10;
f = f(ind);
CS = CS(ind);
Pxy = Pxy(ind);
phase_u1 = unwrap(angle(CS)); %result is in radians
phase_u2 = unwrap(angle(Pxy)); %result is in radians
xlabel('Freq (Hz)');
ylabel('Angle (Rad)');
function [Txy,Cxy,Pxy,f] = tfe_func(x,y,nfft,Fs,window,noverlap)
%TFE Transfer Function Estimate (Txy) / Coherence(Cxy) / Cross spectrum (Pxy)
% applicable to MISO signals (multi input , single output).
% x must be a vector column oriented (m x 1)
% y must be a matrix of dimensions (m x n) (n channels)
% [Txy,Cxy,Pxy,f] = tfe_func(x,y,nfft,Fs,window,noverlap) estimates the transfer function of the
% system with input X and output Y using Welch's averaged periodogram
% method. X and Y are divided into overlapping sections, each of which
% is detrended, then windowed by the WINDOW parameter, then zero-padded
% to length NFFT. The magnitude squared of the length NFFT DFTs of the
% sections of X are averaged to form Pxx, the Power Spectral Density of X.
% The products of the length NFFT DFTs of the sections of X and Y are
% averaged to form Pxy, the Cross Spectral Density of X and Y. Txy
% is the quotient of Pxy and Pxx; it has length NFFT/2+1 for NFFT even,
% (NFFT+1)/2 for NFFT odd, or NFFT if X or Y is complex. If you specify
% a scalar for WINDOW, a Hanning window of that length is used. Fs is
% the sampling frequency which does not effect the transfer function
% estimate but is used for scaling of plots.
% The function returns a vector of freq-
% uencies the same size as Txy at which the transfer function is
% estimated, and overlaps the sections of X and Y by NOVERLAP samples.
% compute PSD and CSD
window = window(:);
n = length(x); % Number of data points
nwind = length(window); % length of window
if n < nwind % zero-pad x , y if length is less than the window length
x = x(:); % Make sure x is a column vector
[p,q] = size(y);
if p~=n & q~=n
error('x an y have incompatible dimensions');
if q==n % y in wrong orientation
y = y';
[p,q] = size(y);
k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;
Pxx = zeros(nfft,1);% Pxx2 = zeros(nfft,1);
Pxy = zeros(nfft,q);% Pxy2 = zeros(nfft,q);
Pyy = zeros(nfft,q);
for i=1:k
xw = window.*x(index);
yw = (window*ones(1,q)).*y(index,:);
index = index + (nwind - noverlap);
Xx = fft(xw,nfft);
Yy = fft(yw,nfft,1);
Xx2 = abs(Xx).^2;
Xy2 = Yy.*(conj(Xx)*ones(1,q));
Yy2 = abs(Yy).^2;
Pxx = Pxx + Xx2;
Pxy = Pxy + Xy2;
Pyy = Pyy + Yy2;
% Select first half
if ~any(any(imag([x y])~=0)), % if x and y are not complex
if rem(nfft,2), % nfft odd
select = [1:(nfft+1)/2];
select = [1:nfft/2+1]; % include DC AND Nyquist
Pxx = Pxx(select);
Pxy = Pxy(select,:);
Pyy = Pyy(select,:);
select = 1:nfft;
% set up output parameters
Txy = Pxy ./ (Pxx*ones(1,q)); % transfer function estimate
Cxy = (abs(Pxy).^2)./((Pxx*ones(1,q)).*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;

More Answers (0)




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!