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

Tags: ftp MATLAB Lambda C

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