用巴特沃斯滤波器进行潮汐滤波分析

作业记录
题目：利用某潮位站一月份的逐时潮汐观测数据，采用巴特沃斯低通滤波器进行潮汐滤波分析，求出低通滤波结果和高通滤波结果。

一、巴特沃斯滤波器及滤波器设计

1. 巴特沃斯滤波器
   巴特沃斯滤波器是一种递归滤波器，它的输出是由输入数据和输出的过去值组成。令ξ=ξ(ω)是频率为ω的正弦和余弦的单调递增有理函数。此单调递增函数生成一个特别有用的具有截止频率ω_c的理想低通递归滤波器的平方增益的近似值：
   
   与非递归滤波器不同的是，巴特沃斯滤波器记录两端没有数据损失，有N个输入数据就能生成N个输出值。但也有个缺陷，就是振铃效应会使滤波输出端的数据失真。因此，像处理非递归滤波器产生的数据缺失问题一样，暂且忽略滤波后末端的输出值。
   巴特沃斯滤波法属于具有时域方程的物理可实现的递归滤波器。因为单个脉冲输入的影响能在未来任意时刻预测到，因此也可将其归类为无限脉冲响应滤波器。高通和带通巴特沃斯滤波器都可以用低通滤波器构建。

2. 滤波器设计
   基于采样间隔 开始指定采样频率ω_s=2πf_s=2π/Δt，这里
   0<ω/ω_s<0.5
   上式中的上限表示归一化奈奎斯特频率ω_N/ω_s。接下来指定在滤波半功率点需要的截止频率ω_C。为了取得最佳效果，滤波器的归一化截止频率ω_c/ω_s应使得滤波器的过渡带与任何采样端区域没有明显的重叠。一旦确定了归一化截止频率，滤波响应的性质完全由滤波器的阶数q 决定。经验表明，参数q应该小于10，最好不大于8。尽管在整个计算中使用双精度，但舍入误差和振铃效应可能会扭曲大值q的滤波器响应，致使滤波器不可用。
   一旦确定好截止频率，巴特沃斯滤波器有两种设计方法。第一种是指定q，使得通带和阻带的衰减能够自动确定。第二种是利用严格单调函数，根据给定频率下所需的衰减计算q。假设我们希望在低通滤波器的阻带有一个截止频率ω_C<ω_a，在频率ω_a处有-D分贝的衰减。利用分贝定义和方程（1），我们发现
   
   其中D是正数，用来计算以分贝为单位的滤波器振幅。如果参数（ω_a,D）已经被正确指定，并且q小于10，那么可以选取最接近的整数取作为滤波器阶数。如果不满足q小于10，那么说明施加的约束条件过于严苛，需要设置新的参数。以上的过程适用于高通滤波器，截止频率在频率ω_a<ω_c衰减为-D时指定q，并将方程（2）中的log⁡〖(ξ_a/ξ_c)〗用log⁡〖(ξ_c/ξ_a)〗替代就能完成。因为log⁡〖(x)〗=-log⁡〖(1/x)〗，忽略log⁡〖(1/x)〗前面的负号，我们能将应用方程（2）简化为高通滤波器。
   获得低通和高通巴特沃斯滤波器的流程图如下图所示。
   
   二、程序实现
   1.代码如下：

clc ;clear all;
%---------原始数据的加载---------%
load('Butterworth.mat', 'data');
d=data';
y=reshape(d,1,[]);

%%%输入信号
fs = 100; %采样频率
Fs=100; %采样频率
Ts=1/Fs; %采样周期 
N=744; %采样个数
t = 1:31;
tt = 0:(1/24):max(t);
tt = tt(1:744);%时间   
figure(1)
subplot(211)
plot(tt,y);%画出原始潮位站一月份的逐小时的海面测量高度图
% axis tight;
title('原始潮位站一月份的逐小时的海面测量高度图','FontSize',16);
x1=xlabel('时间','FontSize',16);
y1=ylabel('海面测量高度/cm','FontSize',16);
set(gca,'xlim',[0,31]);
%set(handles,'xtick',1:1:31) ;
set(gca,'xticklabel',{'1月1日','1月5日','1月10日','1月15日','1月20日','1月25日','1月30日'});
hold on;

%----Butterworth低通滤波器----%
%%%得到原始信号的频谱函数%%%
fft_y=fft(y);%用fft函数作频谱分析，得到的是0~fk内的频谱
%fftshift_y=fftshift(fft_y);%而用fftshift函数得到-fk/2~fk/2内的频谱
f=linspace(-Fs/2,Fs/2,N);
%ifft_y=ifft(ifftshift(fftshift_h));
%plot(f,abs(fft_h));
subplot(212)
%plot(f,abs(fftshift_y));
plot(f,abs(fft_y));
xlabel('频率 (Hz)','FontSize',16); ylabel('幅值','FontSize',16) 
title('原始信号的频谱函数','FontSize',16);
grid on
hold on

%%%用buttord函数求出阶数以及Wn%%%
%butter函数是求Butterworth数字滤波器的系数，在求出系数后对信号进行滤波时用filter函数
wp=0.1*pi;  %通带截止频率
ws=0.4*pi;  %阻带截止频率
Rp=5;   %通带衰减率
Rs=50;  %阻带衰减率
Wp1=2/Ts*tan(wp/2);        %计算角频率,将模拟指标转换成数字指标
Ws1=2/Ts*tan(ws/2); 

%%%滤波器的性能分析图%%%
[n,Wn]=buttord(Wp1,Ws1,Rp,Rs,'s');  %选择滤波器的最小阶数，求出3db的截止频率
[fft_z,P,K]=buttap(n);                  %创建butterworth模拟滤波器，零点，极点，增益
[Bap,Aap]=zp2tf(fft_z,P,K);
[b,a]=lp2lp(Bap,Aap,Wn);   
[bd,ad]=bilinear(b,a,Fs);      %用双线性变换法实现模拟滤波器到数字滤波器的转换
[q,w]=freqz(bd,ad);            %freqz(b,a,计算点数,采样速率) 
figure(2)
subplot(111);
plot(w*fs/(2*pi),abs(q))
grid on;
xlabel('频率 (Hz)','FontSize',16); ylabel('幅值','FontSize',16) 
title('低通滤波器的性能分析图','FontSize',16); 

%%%除去高频信号
figure(3)
N=length(y);
f=fs*(0:fix(N/2-1))/N;%归一化
fft_y=fft(y);
z=filter(bd,ad,y);
subplot(211);
plot(tt,y);
title('(a)原始信号的波形','FontSize',12);
set(gca,'xlim',[0,31]);
set(gca,'xticklabel',{'1月1日','1月5日','1月10日','1月15日','1月20日','1月25日','1月30日'});

subplot(212);
plot(tt,z);
set(gca,'xlim',[0,31]);
title('(b)低通滤波后信号的波形','FontSize',12); 
set(gca,'xticklabel',{'1月1日','1月5日','1月10日','1月15日','1月20日','1月25日','1月30日'});
suptitle('低通滤波前后波形分析');

figure(4)
%sound(z,fs);
subplot(211);
plot(f,abs(fft_y(1:fix(N/2))));
title('(a)原始信号的频谱','FontSize',12);
xlabel('Hz');
fft_z=fft(z);
subplot(212);
plot(f,abs(fft_z(1:fix(N/2)))); 
title('(b)滤波后的信号频谱','FontSize',12);
xlabel('Hz','FontSize',12);
suptitle('低通滤波前后频率对比');

%----用Butterworth低通滤波器构建高通滤波器----%
G = fft_y-fft_z; %用原始信号的频率减去低频信号的频率
g = ifft(G);
figure(5)
subplot(211);
plot(tt,y);
set(gca,'xlim',[0,31]);
title('原始信号的波形','FontSize',12);
set(gca,'xticklabel',{'1月1日','1月5日','1月10日','1月15日','1月20日','1月25日','1月30日'});
subplot(212);
plot(tt,g);
set(gca,'xlim',[0,31]);
title('高通滤波后信号的波形','FontSize',12); 
set(gca,'xticklabel',{'1月1日','1月5日','1月10日','1月15日','1月20日','1月25日','1月30日'});
suptitle('高通滤波后的波形分析');

figure(6)
%sound(z,fs);
subplot(211);
plot(f,abs(fft_y(1:fix(N/2))));
title('原始信号的频谱','FontSize',12);
xlabel('Hz');
subplot(212);
plot(f,abs(G(1:fix(N/2)))); 
set(gca,'ylim',[0,150000]);
title('滤波后的信号频谱','FontSize',12);
xlabel('Hz','FontSize',12);
suptitle('高通滤波前后频率对比');

2. 程序结果如下：
   （1）原始潮位站一月份的逐小时海面测量高度图及其频谱函数如下：
   
   2）设计的巴特沃斯低通滤波器性能分析图如下：
   
   （3）用低通巴特沃斯滤波器对原始信号进行滤波的前后波形对比如下：
   
   （4）用低通巴特沃斯滤波器对原始信号进行滤波的前后频谱对比如下：
   
   （5）用低通巴特沃斯滤波器来获得高通巴特沃斯滤波器，并对原始信号进行高通滤波的前后波形对比如下：
   
   （6）用高通巴特沃斯滤波器对原始信号进行滤波的前后频谱对比如下：
   
   【为什么作业这么多？不知道。。】