matlab滤波器fdatool,各种类型滤波器设计(fdatool，原理，matlab代码)
數據處理
對于一組數據，只有時間戳和加速度，怎么樣進行傅立葉變換分析？ 參考信號處理內容，首先模擬一組數據進行分析。
以下數據兩個頻率為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)=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?
對應代碼為:
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?kH(z)=\sum_{k=0}^{K} b_{k} z^{-k}H(z)=∑k=0K?bk?z?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)=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代表截斷頻率，代碼如下：
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)=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)?
由于存在：
λ=vttbs\lambda=v t_{t b s}λ=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;
標簽：濾波器,filter,width,matlab,fh,trans,data,fdatool
來源： https://blog.csdn.net/chenshiming1995/article/details/104802212
總結
以上是生活随笔為你收集整理的matlab滤波器fdatool,各种类型滤波器设计(fdatool，原理，matlab代码)的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Vue + Spring Boot 项目
- 下一篇： matlab用蒙特卡洛法计算圆周率