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

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)

[Mathieu NOE]

hello

please attach your data files with the paperclip button

all the best

[Faisal]

Actually I have changed the code a bit. I am posting it again in the comment:

sig1 = readtable("total_1770_1776.txt");
sig2 = readtable("HAK_PHA_2017_11_26_1770_1776.txt");
s1 = table2array(sig1);
s2 = table2array(sig2);
Samples = length(s1);
Ts = mean(diff(s1(:,1)));                   % Sampling Timeperiod
Fs = 1/Ts;                                  % Sampling frequency
f = Fs*(0:(Samples/2))/Samples;
F1 = fft(s1(:,2));                          % Fast Fourier Transform of s1
F2 = fft(s2(:,2));                          % Fast Fourier Transform of s2
Corrs1s2 = (ifft(F1.*conj(F2)));            % Cross-Correlation
CS = fft(Corrs1s2);                         % Cross-spectrum
amp = abs(CS);
amp1 = amp(1:Samples/2+1);

Warning: Integer operands are required for colon operator when used as index.

A2 = angle(CS);
A1 = A2(1:Samples/2+1);

Warning: Integer operands are required for colon operator when used as index.

A = unwrap(A1);
plot(f,A, "b-")
xlabel("Frequency (Hz)")
ylabel("Phase (v)")
title("Window 1770 to 1776")
xlim(0:1)

[Faisal]

I am not sure whether the code is ok to find the phase values correctly at the desired frequencies or not.

[Walter Roberson]

Note: to get rid of the warnings, floor(Samples/2) -- you are encountering the case where Samples is odd, so Samples/2 is giving "something plus 1/2" as a result.

[Faisal]

But all I want are the correct values of unwrapped phase of the cross-spectrum at the frquencies from 0.1 to 0.9 with the spacing of 0.1. Can you please tell me whether my code is ok for that or can you suggest any other thing to get the correct results?

[Walter Roberson]

Okay, live with the warnings if you want. It doesn't matter to me.

Answer 1

[Mathieu NOE]

hello

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 :

clc
clearvars

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

plot(f,phase_u1,'b',f,phase_u2,'rd')
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(nwind)=0;
    y(nwind)=0;  
    n=nwind;
end
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');
end

if q==n         % y in wrong orientation
    y = y';
end
[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;

end

% 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];
    else
        select = [1:nfft/2+1];   % include DC AND Nyquist
    end
    Pxx = Pxx(select);
    Pxy = Pxy(select,:);
    Pyy = Pyy(select,:);

else
    select = 1:nfft;
end

% 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;

end

[Faisal]

Thank you very much @Mathieu NOE

[Mathieu NOE]

My pleasure !