# 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=0Kbkz−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πfcTs(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;