How to apply FIR1 Bandpass filter ?

Question

majid 2021-6-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/866765-how-to-apply-fir1-bandpass-filter

评论： majid 2021-6-28

Hi,

I have noisy chirp signal which I want to apply FIR1 Bandpass filter with hanning window to final mixed signal in order to get intermediate frequency signal which I show below and will be derive after applying filtering.

and also may I ask you to explain how we should choose the bandpass filter frequency? frequencies that I used send me errors! and I read help but I cant understande where is my problem in code.

Here is my code:

%%
set(0,'defaultlinelinewidth',2)
clear all
clc
%chirp signal generation
fs=10000;
Ts=1/fs;
t = 0:1/fs:20.5e-3;
f0=20;    f1=320e1;   t1=20.5e-3;% lowered frequencies to show sampling problem
a=(f1-f0)/(t(end)-t(1));  
f=f0+a.*t;
y1 = chirp(t,f0,t1,f1,'linear');
Y = [y1 y1(end:-1:1)]; % y1(end:-1:1) is the time mirrored version of y1, without need to call chirp function again
t1=0:1/fs:41.001e-3;
% figure; plot(t1,Y); title('Bidirectional chirp signal ');
% xlabel('Time(s)'); ylabel('Amplitude');
N=16;
% YY1=repmat(Y,1,N);
% t_YY1=0:1/fs:656.021e-3 ;
% plot(t_YY1,YY1)
% title('Bidirectional chirp signal ');
% xlabel('Time(s)'); ylabel('Amplitude');
clear t_YY1 Y y1 YY1 t1 t
%--------------------------------------------------------------------
%% Delayed signal
R=60;
delta_t=2*R/(3e8);
t_Recieved =delta_t:1/fs:20.5e-3+delta_t;
t_target_Recieved=t_Recieved(end);
y_Recieved =chirp(t_Recieved,f0,t_target_Recieved,f1,'linear');
Y_Recieved = [y_Recieved y_Recieved(end:-1:1)]; 
t_Y_Recieved=delta_t:1/fs:41.001e-3+delta_t;
% figure; plot(t_Y_Recieved.*1000,Y_Recieved); 
YY_Recieved=repmat(Y_Recieved,1,N);
t_YY_Recieved=delta_t:1/fs:659.1e-3+delta_t ;
% plot(t_YY_Recieved,YY_Recieved)
% title('Bidirectional chirp signal ');
% xlabel('Time(s)'); ylabel('Amplitude');
clear delta_t t_Recieved t_target_Recieved t_transmit t_Y_Recieved
clear y_Recieved Y_Recieved
%% adding noise
%speckle noise
M=ones(size(YY_Recieved)); %coherence parameter
x=rand(size(YY_Recieved));
i_speckle=icdf('nbin',x,M,M./(YY_Recieved+M));
clear M x 
% figure;
% plot(t_YY_Recieved.*1000,i_speckle);title(['speckle noise']);
% xlabel('Time (ms)');
% axis tight
% figure
% sig=sqrt(YY_Recieved.^2+i_speckle.^2);
% plot(t_YY_Recieved.*1000,sig);title(['signal + speckle noise']);
% axis tight
%--------------------------------------------------------------------------
%thermal noise
qe=1.602*10e-19; %electron charge
kb=1.3806504*10^(-23);
R_fb=1e6; %feedback resistance (MOhm)
To=300;%absolute temperature
B_det=1.5/(2*Ts);% IF detection band noise bandwidth
quantum_eff=0.2; %detector quantum efficiency
NEPh_dk_pc=1E+6;%Noise equivalent photon rate at the photocathode per sub-pixel from dark current
N_pix_x=14;%number of sub-pixel binned in the x direction to form a super pixel
N_pix_y=14;%number of sub-pixel binned in the y direction to form a super pixel
i_thermal=sqrt(2*qe*kb*To*B_det*N_pix_x*N_pix_y/R_fb)*randn(size(YY_Recieved));%TIA thermal noise
% figure;
% plot(t_YY_Recieved.*1000,i_thermal);
% xlabel('Time (ms)');title(['thermal noise']);
% axis tight
% figure
% sig=sqrt(i_thermal.^2+YY_Recieved.^2);
% plot(t_YY_Recieved.*1000,sig);title(['signal + thermal noise']);
% axis tight
clear kb R_fb To B_det sig
%--------------------------------------------------------------------------
%shot noise modeled by white noise
SNR=20;
i_shot=awgn(YY_Recieved,SNR);% white noise adding
% figure;
% plot(t_YY_Recieved.*1000,abs(i_shot));title(['shot noise modeled by white gaussian noise']);
% xlabel('Time (ms)');
% axis tight
% figure
% sig=sqrt(Y_Recieved.^2+i_shot.^2);
% plot(t_YY_Recieved.*1000,sig);title(['signal + gaussian noise']);
% axis tight
%--------------------------------------------------------------------------
noisy_sig1=i_thermal.^2+i_shot.^2+i_speckle.^2;
clear SNR sig i_thermal i_shot i_speckle x_8
%--------------------------------------------------------------------------
% backreflection , backscattering , atmospheric effects
P_l_pk=7.4; %laser peak power (watt)
FF=0.49; %Reciever FPA fill factor(%)
t_t=0.7; %transmitter transmitance
t_r=0.338; %reciever transmitance
D_r=0.254; %Reciever diameter (m)
rho_tgt=0.1; %target reflectance
kapa_tgt_gauss=1; %gaussian beam target return geometrical form factor
alpha=0.0002; %atmospheric extinction coefficient
rho_bg=0.3; %background reflectance
kapa_ls_gauss=1;%gaussian spot geometrical form factor
beta_vol_bs=9.95e-6;
Hs=267; %solar spectral irradiance (watt/(m^2 * um))
N_sky=3; %sky spectral irradiance (watt/(m^2 * sr*um))
delta_lambda=6.2e-9; %filter bandwidth
delta=12e-5; %reciever sub pixel size
Reciever_focal=1.016; %(m)
A_ifov=(R*delta/Reciever_focal)^2;
D_obst=0.5; %reciver central obstruction diameter
A_r=pi*D_r^2/4;%reciever aperture area
A_obst=pi*D_obst^2/4;%reciever central obstruction area
omega_rcv=(A_r-A_obst)/R^2;%solid angle subtended by reciever aperture at the target
Pb=abs((rho_bg*Hs/pi*exp(-alpha*R)+N_sky)*delta_lambda*t_r*FF*A_ifov*omega_rcv);
%power recieved after narrowband filter
%--------------------------------------------------------------------------
mt=0.5;%transmiter modulation index
mr=0.5;%reciever modulation index
em=1;%mixing efficiency
G=20;%electron gain of tube
R_max=0.243;%peak responsivity of tube(A/Watt)
gamma_mod_bg=sqrt((1-em*mr/2)^2+((em)^2*(mr^2)^2)/2);%
i_b=G*R_max*gamma_mod_bg*Pb;%%current of solar background
i_dk=G*qe*quantum_eff*gamma_mod_bg*NEPh_dk_pc*N_pix_x*N_pix_y; %dark current
i=(i_b+i_dk)*10^8;
i_back=poissrnd(i*ones(size(YY_Recieved)));
clear P_l_pk FF t_t t_r D_r rho_tgt kapa_tgt_gauss alpha rho_bg kapa_ls_gauss beta_vol_bs
clear Hs N_sky di_thermalelta_lambda delta Reciever_focal A_ifov D_obst A_r A_obst omega_rcv Pb
clear mt mr em G R_max gamma_mod_bg i_b i_dk i quantum_eff NEPh_dk_pc N_pix_x N_pix_y qe delta_lambda
% figure;
% plot(t_YY_Recieved.*1000,i_back);
% xlabel('Time (ms)');title(['backscattering noise']);
% axis tight
% 
% figure
% sig=sqrt(YY_Recieved.^2+i_back.^2);
% plot(t_YY_Recieved.*1000,sig);title(['signal + atmospheric effects']);
% axis tight
%--------------------------------------------------------------------------
% noisy signal
noisy_signal=sqrt(YY_Recieved.^2+i_back.^2+noisy_sig1);
% figure;
% plot(t_YY_Recieved.*1000,noisy_signal);title(['noisy signal']);
% xlabel('Time (ms)');
% axis tight
clear  i_back YY_Recieved noisy_sig1
%% LO signal
delta_t_LO=193e-9;
t_LO =delta_t_LO:1/fs:20.5e-3+delta_t_LO;
t_target2=t_LO(end);
y_LO =chirp(t_LO,f0,t_target2,f1,'linear');
Y_LO = [y_LO y_LO(end:-1:1)]; 
t_Y_LO=delta_t_LO:1/fs:41.001e-3+delta_t_LO;
% figure; plot(t_Y_LO.*1000,Y_LO); 
YY_LO=repmat(Y_LO,1,N);
t_YY_LO=delta_t_LO:1/fs:659.1e-3+delta_t_LO ;
% plot(t_YY_LO,YY_LO)
% title('Bidirectional chirp signal ');
% xlabel('Time(s)'); ylabel('Amplitude');
clear delta_t_LO t_LO  t_target2 y_LO t_Y_LO t_YY_Recieved Y_LO
%% signal mixing
mixed_signal=YY_LO.*noisy_signal;
plot(t_YY_LO,mixed_signal);
axis tight
clear YY_LO noisy_signal