数据处理与滤波器性质分析

1 问题描述

对于一组数据，**只有时间戳和加速度，怎么样进行傅立叶变换分析？**参考信号处理内容，首先模拟一组数据进行分析。
以下数据两个频率为1Hz与100Hz，经过采样和傅立叶变化之后，捕捉到信号对应的频率为1Hz与100Hz(还有其他信号)。

clear all;
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);        % 滤波后的信号

信号时域频域的关系如下：



因此经常设计的滤波器一般有如下形式：

对应代码为：

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,samplerate);

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

2 滤波器
2.1 FIR滤波器
特点如下：

对于上述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），如下：

总体来说，是下面三个流程：

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

设计过程，可以参考下面：

2.2 IIR 无限滤波器


3 滤波器设计（参考轨检）
容易想到的是，在这里做的数据的卷积处理，放在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;

3.1 模拟滤波器与数字滤波
模拟滤波器如下所示：

二阶低通滤波器代码如下，该滤波器是从模拟滤波器转换而来。

% 二阶低通滤波器
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;

3.2 模拟滤波器与数字滤波器的推导1
4 小波变换
当信号随着时间发生变化时，可能信号的频率随着时间在不断增大，如何观测信号中的频率？其中低频的层粉需要较长的时间测量。
大概得到如下的结果：

https://blog.csdn.net/luoshi006/article/details/51459884 ↩︎