Determine filter function from input and output signals
조회 수: 6 (최근 30일)
이전 댓글 표시
I have an known input signal that I am putting through a "black box" resulting in an experimentally measured output signal. I want to model the "black box" by deducing its filter function from the input and output signal. How can I do that?
The code below is a test code to develop the filter deduction capabilities. First I generate a square wave. Then I create the output signal using a first order butterworth filter. This is to simulate the "black box". Then I use invfreq to try to deduce the filter form. However, it does not recreate the orignal butterworth filter.
clear
close all
%% Generate Signal
%sampling frequency
f_sample = 200e6; %samples/s
%sampling period
tau_sample = 1/f_sample; %s
%time vector
t = 0:tau_sample:500e-6; %s
%vector length
L = length(t);
%frequencies
f = (f_sample*(0:(L/2))/L)';
%wave period
tau_wave = 100e-6; %s
%input square wave
vin = square(t*2*pi/(tau_wave));
%% Filter Input signal through "Black Box" to create output signal
%nyquist sampling frequency
Ny = f_sample/2;
%order
n=1;
%cutoff frequency
LP = 14e3;
Wn = LP/Ny;
[b1,a1] = butter(n, Wn,'low');
H1=dfilt.df2t(b1,a1);
HGfilter=H1;
%frequency response of filter
[h,w] = freqz(b1,a1,length(f),f_sample);
%output signal after filter
vout = filter(HGfilter,vin);
%% Fourier Transform
% raw FFT
Vin = fft(vin);
%normalize FFT coefficients
P2 = abs(Vin/L);
Pin = P2(1:L/2+1);
Pin(2:end-1) = 2*Pin(2:end-1);
% raw FFT
Vout = fft(vout);
%normalize FFT coefficients
P2 = abs(Vout/L);
Pout = P2(1:L/2+1);
Pout(2:end-1) = 2*Pout(2:end-1);
%% Deduced filter
%filter transmission
T = abs(Pout./Pin);
%calc transfer function coefficients
[b,a] = invfreqz(T,f,2,3);
%transfer funciton
hz = tf(b,a);
%frequency response
[h_calc, w_calc] = freqz(b,a,length(f),f_sample);
%model filtered signal
[v_calc,t_calc] = lsim(hz,vin,t);
% raw FFT
Vcalc = fft(v_calc);
%normalize FFT coefficients
P2 = abs(Vcalc/L);
Pcalc = P2(1:L/2+1);
Pcalc(2:end-1) = 2*Pcalc(2:end-1);
%% Plot
figure(101)
set(gca, 'FontSize', 18);
plot(t, vin, 'LineWidth', 2)
hold on
plot(t, vout, 'LineWidth', 2);
plot(t_calc,v_calc, 'LineWidth', 2);
hold off
%axis([2.42 2.52 -0.7 10]);
xlabel('Time (\mus)','FontSize', 20, 'FontWeight', 'bold');
ylabel('Signal (a.u.)', 'FontSize', 20, 'FontWeight', 'bold');
legend('Input', 'Measured', 'Filtered');
figure(102)
set(gca, 'FontSize', 18);
plot(f, Pin, 'LineWidth', 2)
hold on
plot(f, Pout, 'LineWidth', 2);
plot(f, Pcalc, 'LineWidth', 2);
hold off
%axis([2.42 2.52 -0.7 10]);
xlabel('Frequency (Hz)','FontSize', 20, 'FontWeight', 'bold');
ylabel('FT (a.u.)', 'FontSize', 20, 'FontWeight', 'bold');
legend('Input', 'Measured','Filtered');
set(gca,'XScale','log')
set(gca,'YScale','log')
%plot filter response
figure(103)
plot(w,h)
hold on
plot(w_calc,h_calc)
hold off
xlabel('Frequency (Hz)','FontSize', 20, 'FontWeight', 'bold');
ylabel('Transfer Function (a.u.)', 'FontSize', 20, 'FontWeight', 'bold');
legend('True Filter','Deduced Filter')
set(gca,'XScale','log')
set(gca,'YScale','log')
댓글 수: 0
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Signal Generation and Preprocessing에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!