# How can I find the Transfer Function having Magnitude(dB), Phase(degree) and Frequency(Hz)?

20 views (last 30 days)
Liang Kar Yan on 10 Dec 2021
Commented: Liang Kar Yan on 14 Dec 2021
I have 3 individual files which are Magnitude(dB), Phase(degrees) and Frequency(Hz) in excel. I need to find the transfer function with these data. After getting the transfer function, I need to plot back the graph (magnitude and phase) from transfer function to compare with my data.
I wish to get something like the picture attached below.
Liang Kar Yan on 14 Dec 2021
Hi @Star Strider, thank you for your help.

Mathieu NOE on 14 Dec 2021
hello
I tried a few options , IIR or FIR filters fit.
As the phase plot shows , there is a quite significant phase roll rate, idicating the presence of a huge delay in the system.
I assumed a sampling rate of Fs = 1000 hz and found out that more or less we can fit either a FIR or a IIr high pass filter in series with almost 200 samples of delay (hudge !!)
maybe there are more powerful tools then the simple invfreqz (not giving any good results here) or the manual fit I am doing here
FIR fit plot :
IIR fit plot :
clc
clearvars
figure(1)
subplot(211),plot(Freq,Magn_dB);
subplot(212),plot(Freq,Phse);
Magn = 10.^(Magn_dB/20) ;
%% high pass filter model (IIR)
Fs = 1000; % ? to be confirmed
N = 8; % filter order
dc_gain = Magn(end); % asymptotic value
[val,ind] = min(abs(Magn_dB - Magn_dB(end) + 5)); % - 5dB (vs dc_gain) cut off frequency index search
fc = Freq(ind); % - 5dB (vs dc_gain) cut off frequency
[b,a] = butter(N,2*fc/Fs,'high');
b = b*dc_gain; % apply dc gain on numerator
[g,p] = dbode(b,a,1/Fs,2*pi*Freq);
% adding delay due to sampling
nd = 200; % delay (samples)
rpd = -360*nd*Freq/Fs;
p = p+rpd; % adding filter phase to samples delay phase
p = mod(p,360);
p = p -180; % polarity correction
figure(1)
subplot(211),plot(Freq,Magn_dB,Freq,20*log10(g));
title('IIR model fit')
ylabel('Modulus (dB)');
subplot(212),plot(Freq,Phse,Freq,p);
xlabel('Frequency (Hz)');
ylabel('Phase (°)');% return
%% high pass filter model (FIR)
Fs = 1000; % ? to be confirmed
N = 20;
dc_gain = Magn(end); % asymptotic value
[val,ind] = min(abs(Magn_dB - Magn_dB(end) + 3)); % - 3dB (vs dc_gain) cut off frequency index search
fc = Freq(ind); % - 3dB (vs dc_gain) cut off frequency
[b,a] = fir1(N,2*fc/Fs,'high');
b = b*dc_gain; % apply dc gain on numerator
[g,p] = dbode(b,a,1/Fs,2*pi*Freq);
% adding delay due to sampling
Fs = 1000; % ? to be confirmed
nd = 200; % delay (samples)
rpd = -360*nd*Freq/Fs;
p = p+rpd; % adding filter phase to samples delay phase
p = mod(p,360);
p = p -180; % polarity correction
figure(2)
subplot(211),plot(Freq,Magn_dB,Freq,20*log10(g));
title('FIR model fit')
ylabel('Modulus (dB)');
subplot(212),plot(Freq,Phse,Freq,p);
xlabel('Frequency (Hz)');
ylabel('Phase (°)');
%% Id with invfreqz (FIR)
h = Magn .* exp(1j*180/pi*(Phse));
nb = 40+nd;
na = 1;
iter = 1000;
[bb,aa] = invfreqz(h,pi*Freq/Fs,nb,na,[],iter); % stable approximation to system
[g,p] = dbode(bb,aa,1/Fs,2*pi*Freq);
p = mod(p,360);
figure(3)
subplot(211),plot(Freq,Magn_dB,Freq,20*log10(g));
subplot(212),plot(Freq,Phse,Freq,p);
Liang Kar Yan on 14 Dec 2021
Thanks for the useful information.
Thank you.