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

20 views (last 30 days)
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.
  5 Comments

Sign in to comment.

Accepted Answer

Mathieu NOE
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
Freq = readmatrix('Frequency.xlsx');
Magn_dB = readmatrix('Magnitude.xlsx');
Phse = readmatrix('Phase.xlsx');
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);
  10 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!