数据处理
对于一组数据，只有时间戳和加速度，怎么样进行傅立叶变换分析？ 参考信号处理内容，首先模拟一组数据进行分析。
以下数据两个频率为1Hz与100Hz，经过采样和傅立叶变化之后，捕捉到信号对应的频率为1Hz与100Hz(还有其他信号)。
close all;
t = 0:0.01:3; % 真实世界时间
f1 = 1; % 频率
f2 = 200;
f3 = 50; % 设定两个复信号
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)); % 信号函数
y = F(t); % 生成信号
% figure;subplot(3,1,1);plot(t , y); % 信号真实图
fs = 1000; % 采样率
dtc = 1/fs; % 采样间隔时间
tc = 0:dtc:4; % 采样时间序列
yc = F(tc); % 采样信号序列
%% 傅立叶变换以及画图
figure;
N = length(yc);
x = (-N/2+1:N/2)/N*fs;
semilogy(x , abs(fftshift(fft(yc))));
我们可以看到，复信号在幅度谱上表现是只有单侧有信号。而实信号在幅度谱上两侧均有信号。
那么如何对数据进行信号处理呢？如何用fdatool设计滤波器？
频域上表现如下：
设计上述高通滤波器，与所有数据进行卷积，完成滤波。得到结果如下：
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); % 滤波后的信号
信号时域频域的关系如下：
因此经常设计的滤波器一般有如下形式：
H(z)=1−0.2z−1+0.8z−20.2+0.5z−1
对应代码为:
clear, close all
%% initialize parameters
% 载波频率
samplerate = 1000; % in Hz 采样率
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
% % 因为相位的变化会带来一定的相位偏移
% phase = angle(H)*2*pi; % get phase in deg.
%
% faxis = samplerate/2*linspace(0,1,N/2); % the axis of frequency
%% 或者下面：
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滤波器
特点如下：
转换函数为：
H(z)=∑k=0Kbkz−k
对于上述fdatool设计的FIR滤波器，a为0，所以只用b进行卷积运算。下面画出了相位谱和幅度谱，下面作为示例。
%% 设计滤波器（FIR）
N = 512;
a = 1;
H = fft(b,N)/fft(a,N); % H矩阵
mag = 20*log10(abs(H)); % get magnitude of spectrum in dB 幅值
phase = angle(H)*2*pi; % get phase in deg.相位
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
滤波器设计离不开这个函数，具有特殊性质的函数sinc（t），如下：
所以设计以下低通滤波器：
b(k)=π(k−L/2)sin[2πfcTs(k−L/2)]
fc代表截断频率，代码如下：
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))))
% 加窗
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
设计过程，可以参考下面：
那么如何利用matlab代码生成滤波器？
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 无限滤波器
%%% 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');
%% 简单如下
%% 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')';
小波变换
当信号随着时间发生变化时，可能信号的频率随着时间在不断增大，如何观测信号中的频率？其中低频的层粉需要较长的时间测量。
大概得到如下的结果：
滤波器设计
容易想到的是，在这里做的数据的卷积处理，放在c语言中肯定是不合理的。那么在轨检模型中是如何完成计算的？怎么样与之同步起来？
下面给出了两个滤波器设计：
% FMIctrl中的滤波器幅频频特性
% ---------- 10 Hz(对于什么？) -------
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;
模拟滤波器与数字滤波
模拟滤波器如下所示：
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)
由于存在：
λ=vttbs
二阶低通滤波器代码如下，该滤波器是从模拟滤波器转换而来。
% 二阶低通滤波器
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 ('二阶数字抗混叠滤波器和补偿滤波器');
semilogx(v1./f2,20*log10(abs(h2)));hold on;