difference in transfer function calculation in matlab and python
Show older comments
i have developed below code in matlab and trying implement samething in python
but getting different outputs from matlab and python.
can you help here to convert this code to python
Matlab code
[Y, freq] = tfestimate(in, out, hanning(NFFT), round(overlap*NFFT), [], sample_rate);
example output
freq:freq= [0;0.0101800266860033]
Python code
f, Pxy = csd(out_signal,in_signal, fs=sample_rate, window=windows.hann(NFFT, sym=True), average='median', nperseg=NFFT, noverlap=int(overlap * NFFT), scaling='spectrum')
example output
freq: [0; 0.0132058240]
3 Comments
Mathieu NOE
on 11 Apr 2025
hello
seems to me you're comparing apples and oranges
csd gives the cross power spectral density , not a transfer function estimate
fyi : this is the code (of an older version) of tfestimate ; in the litterature this is called the H1 estimator and the Txy transfer function uses Pxx and Pxy :
for zero extra fee you get also the coherence as you have also Pyy :
%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = mytfe_and_coh(x,y,nfft,Fs,window,noverlap)
% Transfer Function and Coherence Estimate
% 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
y = y(:); % Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;
Pxx = zeros(nfft,1);
Pyy = zeros(nfft,1);
Pxy = zeros(nfft,1);
for i=1:k
xw = window.*x(index);
yw = window.*y(index);
index = index + (nwind - noverlap);
Xx = fft(xw,nfft);
Yy = fft(yw,nfft);
Xx2 = abs(Xx).^2;
Yy2 = abs(Yy).^2;
Xy2 = Yy.*conj(Xx);
Pxx = Pxx + Xx2;
Pyy = Pyy + Yy2;
Pxy = Pxy + Xy2;
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);
Pyy = Pyy(select);
Pxy = Pxy(select);
else
select = 1:nfft;
end
Txy = Pxy ./ Pxx; % transfer function estimate
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;
end
siva prasad
on 14 Apr 2025
Mathieu NOE
on 14 Apr 2025
hello again
as I said before , the python code you mention is not the transfer estimate
Answers (0)
Categories
Find more on Call Python from MATLAB in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!