Various types of filter design (fdatool, principle, matlab code)

data processing

For a group of data, only time stamp and acceleration, how to do Fourier transform analysis? Firstly, a set of data is simulated for analysis.
The following data have two frequencies of 1Hz and 100Hz. After sampling and Fourier transform, the frequency corresponding to the signal captured is 1Hz and 100Hz (and other signals).

	close all;
	
	t = 0:0.01:3;               % Real world time
	f1 = 1;                      % frequency
	f2 = 200;           
	f3 = 50;                    % Set two complex signals
	f4 = -60;
	
	F = @(t)(sin(2*pi*f1*t) + sin(2*pi*f2*t)+ exp(j*2*pi*f3*t) + exp(j*2*pi*f4*t));    % Signal function
	y = F(t);                   % Generating signal
	% figure;subplot(3,1,1);plot(t , y);         % True signal diagram
	fs = 1000;                  % sampling rate
	dtc = 1/fs;                 % Sampling interval
	tc = 0:dtc:4;               % Sampling time series
	yc = F(tc);          % Sampling signal sequence
	%% Fourier transform and drawing
	figure;
	N = length(yc);
	x = (-N/2+1:N/2)/N*fs;
	semilogy(x , abs(fftshift(fft(yc))));

We can see that only one side of complex signal has signal in amplitude spectrum. The real signal has signals on both sides of the amplitude spectrum.
So how to process the data? How to design filter with fdatool?


The frequency domain is as follows:

The high pass filter is designed and convoluted with all data to complete the filtering. The results are as follows:

	Fs = 1000;  % Sampling Frequency
	Fstop = 50;              % Stopband Frequency
	Fpass = 100;             % Passband Frequency
	Dstop = 0.0001;          % Stopband Attenuation
	Dpass = 0.057501127785;  % Passband Ripple
	dens  = 20;              % Density Factor
	
	% Calculate the order from the parameters using FIRPMORD.
	[N, Fo, Ao, W] = firpmord([Fstop, Fpass]/(Fs/2), [0 1], [Dstop, Dpass]);
	
	% Calculate the coefficients using the FIRPM function.
	b  = firpm(N, Fo, Ao, W, {dens});
	
	Hd = dfilt.dffir(b);
	
	yf = conv( b , yc);		% Filtered signal

The relationship between signal time domain and frequency domain is as follows:

Therefore, filters often designed are generally in the following forms:
H(z)=0.2+0.5z−11−0.2z−1+0.8z−2H(z)=\frac{0.2+0.5 z^{-1}}{1-0.2 z^{-1}+0.8 z^{-2}}H(z)=1−0.2z−1+0.8z−20.2+0.5z−1​

The corresponding code is:

clear, close all

%% initialize parameters
% carrier frequency
samplerate = 1000;   % in Hz sampling rate
N = 512;             % number of points, must be even, better be power of 2

%% define a and b coeffients of H (transfer function)
a = [1 -0.2 0.8];   % denominator terms
b = [0.2 0.5];      % numerator terms

%% option 1:compute the spectrum of H using fft
% H = fft(b,N)./fft(a,N);    % compute H(f)
% 
% mag = 20*log10(abs(H));    % get magnitude of spectrum in dB
% % Because the change of phase will lead to a certain phase shift
% phase = angle(H)*2*pi;     % get phase in deg.
% 
% faxis = samplerate/2*linspace(0,1,N/2);  % the axis of frequency
%% Or as follows:
N = 512;
[h1 , ftp] = freqz(b,1,N,fs);

mag = 20*log10(abs(h1));    % get magnitude of spectrum in dB
phase = angle(h1)/pi*180;     % get phase in deg.

figure,
subplot(2,1,1),plot(ftp,mag)
xlabel('Frequency (Hz)'),ylabel('Magnitude (dB)')
grid on

subplot(2,1,2),plot(ftp,phase,'r')
xlabel('Frequency (Hz)'),ylabel('Phase (deg.)')
grid on


FIR filter

The features are as follows:

The conversion function is:
H(z)=∑k=0Kbkz−kH(z)=\sum_{k=0}^{K} b_{k} z^{-k}H(z)=∑k=0K​bk​z−k
For the FIR filter designed by fdatool, a is 0, so only b is used for convolution. The phase spectrum and amplitude spectrum are drawn below, as an example below.

		%% Design filter( FIR)
		N = 512;
		a = 1;
		H = fft(b,N)/fft(a,N);   % H matrix
		mag = 20*log10(abs(H));    % get magnitude of spectrum in dB amplitude
		phase = angle(H)*2*pi;     % get phase in deg.phase
		faxis = samplerate/2*linspace(0,1,N/2);  % the axis of frequency
		%% plot the spectrum of H
		figure,
		subplot(2,1,1),plot(faxis,mag(1:N/2))
		xlabel('Frequency (Hz)'),ylabel('Magnitude (dB)')
		grid on
		subplot(2,1,2),plot(faxis,phase(1:N/2),'r')
		xlabel('Frequency (Hz)'),ylabel('Phase (deg.)')
		grid on

This function is indispensable for filter design. sinc (t), which has special properties, is as follows:


Therefore, the following low-pass filters are designed:
b(k)=sin⁡[2πfcTs(k−L/2)]π(k−L/2)b(k)=\frac{\sin \left[2 \pi f_{c} T_{s}(k-L / 2)\right]}{\pi(k-L / 2)}b(k)=π(k−L/2)sin[2πfc​Ts​(k−L/2)]​
fc stands for truncation frequency, the code is as follows:

		L = 57;
		fs = 1000;
		f2 = 100;
		for k = 1:L
		    b(k) = sin(2*pi*f2*dtc*(k - L/2))/(pi*(k-L/2));
		end
		figure;
		N = length(b);
		x = (-N/2+1:N/2)/N*fs;
		semilogy( x,abs(fftshift(fft(b))))

		% Windowing
		faxis = fs/2*linspace(0,1,N/2);
		HW = fft(b.*hamming( length(b) )',N);
		mag = 20*log10(abs(HW));
		figure
		plot(faxis,mag(1:N/2))
		xlabel('Frequency (Hz)'),ylabel('Magnitude (dB)')
		grid on

For the design process, please refer to the following:

So how to use matlab code to generate filter?

fl=75;   % low-cutoff frequency
		fh=165;  % high-cutoff frequency
		trans_width=20; % in Hz. It is a half of transition band. if data length is not long enough, increase trans_width
		
		rp=1;  % in dB
		rs=40; % in dB
		
		%%% lowpass filter
		[data_3sFIR,forder] = filter_3sFIR(data,[fl-trans_width fl+trans_width],[1 0],[0.1 0.001],samplerate); 
		
		%%% bandpass filter
		[data_3sFIR,forder] = filter_3sFIR(data,[fl-trans_width fl+trans_width fh-trans_width fh+trans_width],[0 1 0],[0.001 0.1 0.001],samplerate); 
		
		%%% highpass filter
		[data_3sFIR,forder] = filter_3sFIR(data,[fh-trans_width fh+trans_width],[0 1],[0.001 0.1],samplerate); 
		
		%%% bandstop filter
		[data_3sFIR,forder] = filter_3sFIR(data,[fl-trans_width fl+trans_width fh-trans_width fh+trans_width],[1 0 1],[0.1 0.001 0.1],samplerate); 

IIR infinite filter

%%% lowpass filter
		[data_3sIIR,forder] = filter_3sIIR(data,fl-trans_width,fl+trans_width,rp,rs,samplerate,'low'); 
		
		%%% bandpass filter
		[data_3sIIR,forder] = filter_3sIIR(data,[fl+trans_width fh-trans_width],[fl-trans_width fh+trans_width],rp,rs,samplerate,'bandpass'); 
		
		%%% highpass filter
		[data_3sIIR,forder] = filter_3sIIR(data,fh+trans_width,fh-trans_width,rp,rs,samplerate,'high'); 
		
		%%% bandstop filter
		[data_3sIIR,forder] = filter_3sIIR(data,[fl-trans_width fh+trans_width],[fl+trans_width fh-trans_width],rp,rs,samplerate,'stop');  
		%% Simple as follows
		%% filter 
		sigfilter1=filter_2sIIR(EEGdata',fh,samplerate,forder,'low')';
		
		sigfilter2=filter_2sIIR(EEGdata',fl,samplerate,forder,'high')';
		
		sigfilter3=filter_2sIIR(EEGdata',[fl fh],samplerate,forder,'bandpass')';

wavelet transform

When the signal changes with time, the frequency of the possible signal increases with time. How to observe the frequency in the signal? The low frequency layer powder needs a long time to measure.
The results are as follows:

filter design

It is easy to think that the convolution processing of data done here is certainly unreasonable in c language. So how to complete the calculation in orbit inspection model? How to synchronize with it?
Two filter designs are given below:

% FMIctrl Amplitude frequency characteristics of filters in
% ---------- 10 Hz(For what?)  -------
fs = 500;
N = 80000;
b10 = [40000 0 0];
a10 = [4010000 -7600000 3610000];
[h10 f10]= freqz(b10,a10,N,'whole',fs);

% 
mag = 20*log10(abs(h10));    % get magnitude of spectrum in dB
phase = angle(h10)/pi*180;     % get phase in deg.
figure,
subplot(2,1,1),semilogx(f10,mag)
xlabel('Frequency (Hz)'),ylabel('Magnitude (dB)')
grid on
subplot(2,1,2),semilogx(f10,phase,'r')
xlabel('Frequency (Hz)'),ylabel('Phase (deg.)')
grid on
suptitle('10Hz');

% ----------20 Hz-----------
coef1 = 40000;coef2= 1800000;
coef3=810000 ;coef4=10340000 ;
b20 = [coef1 0 0];
a20 = [coef4 -coef2 coef3];
figure();
[h20 f20]= freqz(b20,a20,N,'whole',fs);
subplot(2,1,1);semilogx(f20,20*log10(abs(h20)));xlabel('Frequency (Hz)'),ylabel('Magnitude (dB)')
subplot(2,1,2);semilogx(f20,angle(h20)*180/pi);xlabel('Frequency (Hz)'),ylabel('Phase (deg.)')
suptitle('20Hz');grid on;

Analog filter and digital filter

The analog filter is as follows:
H(s)=B(s)A(s)=b(1)sn+b(2)sn−1+⋯+b(n+1)a(1)sm+a(2)sm−1+⋯+a(m+1) H(s)=\frac{B(s)}{A(s)}=\frac{b(1) s^{n}+b(2) s^{n-1}+\dots+b(n+1)}{a(1) s^{m}+a(2) s^{m-1}+\dots+a(m+1)} H(s)=A(s)B(s)​=a(1)sm+a(2)sm−1+⋯+a(m+1)b(1)sn+b(2)sn−1+⋯+b(n+1)​
Due to the presence:
λ=vttbs \lambda=v t_{t b s} λ=vttbs​
The second-order low-pass filter code is as follows, which is converted from the analog filter.

% Second order low pass filter
w2 = (10^5)/(2^14);
v1=  15/3.6;
t1= 0.25/v1;
w2t1 = w2*t1;
b2 = [(w2t1)^2 0 0];
a2 = [1+w2t1+(w2t1)^2  ,- (2 + w2t1)  ,1];
[h2 f2] = freqz(b2,a2,800000,500);
figure;	suptitle ('Second order digital anti aliasing filter and compensation filter');
semilogx(v1./f2,20*log10(abs(h2)));hold on;
Published 4 original articles, won praise 0, visited 170
Private letter follow

Tags: ftp MATLAB Lambda C

Posted on Wed, 11 Mar 2020 06:50:25 -0400 by gerkintrigg