MATLAB仿真分别实现一阶RC低通和高通滤波器，输入信号为正弦信号或者方波信号。
注意截止频率为f = 1/(2*pi*R*C)
低通滤波器下所示：
%功能：一阶RC低通滤波器仿真 %说明： %1、分析了一阶RC滤波器的幅值衰减特性和相移特性 %2、分析了一阶RC滤波器的频域特性 %3、使用lsim对系统进行仿真 %4、使用FFT对原始输入信号和滤波器输出信号进行分析 %传递函数：sys=1/(1+sRC) %========================================================================== clc clear close all A=1; % A 幅度值 fs=5000000; % fs 采样率//最后转化为频域的最远距离最大频率//采样的频率 F=10000; % F 频率,小于采样率的一半（奈奎斯特）//信号频率 N=5000; % N 采样个数//取了多少个点// %采样频率与采样个数没有什么必然关系，采样个数乘上采样频率的倒数（时间间隔）即为总的时间长度 dt=1/fs; %时间间隔//决定了时间点与点之间的距离，当然随便指定一个数也行，只不过使信号的时域上的长度不同 t=0:dt:(N-1)*dt; %时间向量 F1=1.0*10^3; F2=5.0*10^3; F3=8.0*10^3; F4=25*10^3; F5=45*10^3; F6=65*10^3; F7=100*10^3; F8=300*10^3; F9=400*10^3; F10=15.9*1000; %------------------------信源产生------------------------------------------- dataSourceType=0; %波形类型 0表示正弦波，1表示三角波，2表示方波 switch dataSourceType case 0 % y=A*sin(2*pi*F1*t); % y=A*sin(2*pi*F2*t); % y=A*sin(2*pi*F3*t); % y=A*sin(2*pi*F4*t); % y=A*sin(2*pi*F5*t); % y=A*sin(2*pi*F6*t) % y=A*sin(2*pi*F7*t); % y=A*sin(2*pi*F8*t); % y=A*sin(2*pi*F9*t); y=A*sin(2*pi*F10*t); case 1 y=A*sawtooth(2*pi*F*t,0.5); %三角波 case 2 y=A*square(2*pi*F*t,50); %方波 otherwise end %--------------------画出原始输入信号的时域与频域图像---------------------- figure(1); subplot(2,1,1); plot(t,y); title('信号的时域波形');%对原始图像进行时域画图 xlabel('时间/s'); ylabel('电压/v'); h=fft(y,5000);%快速傅里叶变换 h_d=abs(fftshift(h));%使频域图像中间为零 f=(-N/2:N/2-1)*fs/N;%将取得时间上的点转化为频率上的点 subplot(2,1,2); plot(f,h_d/5000);%画原始图像频域上图 %axis([-4*10^4 4*10^4 0 1]); title('信号的频域波形'); xlabel('频率/hz'); ylabel('电压/v'); %---------------------------滤波器部分------------------------------------ %-----------------------滤波器本身特性（波特图）---------------------------- r=10.2*10^3;%电阻阻值（ome） c=((1200*5600)/(1200+5600))*10^(-12);%电容（f） w=2*pi*f;%角频率 Para=r*c*1i; for n=1:length(f) S(n)=abs(1/(1+Para*w(n)));%滤波器对于不同频率幅值衰减系数 P(n)=angle(1/(1+Para*w(n)))*180/pi;%滤波器对于不同频率相移系数 end figure(2); subplot(2,1,1); plot(f,S,'r');%幅值曲线 title('幅值衰减特性'); xlabel('频率/hz'); ylabel('增益倍数'); subplot(2,1,2); plot(f,P,'blue');%相位曲线 title('相位特性'); xlabel('频率/hz'); ylabel('相位'); Func=tf(1,[r*c,1]);%系统的传递函数 figure(3); bode(Func);%系统的波特图 title('幅频特性'); %--------------------信号通过滤波器---------------------------------- [yout,tout] = lsim(Func,y,t);%滤波后信号图像 %-----------------------时域图像------------------------------------ figure(4); subplot(2,1,1); plot(t,y); title('原始信号'); xlabel('时间/s'); ylabel('电压/v'); subplot(2,1,2); plot(tout,yout); title('滤波后的时域波形'); xlabel('时间/s'); ylabel('电压/v'); %------------------------------频域图像--------------------------------- q=fft(yout); q_d=abs(fftshift(q)); figure(5); subplot(2,1,1); plot(f,q_d); title('滤波后的频谱'); xlabel('频率/hz'); ylabel('电压/v'); subplot(2,1,2); plot(f,h_d); title('输入信号的频谱'); xlabel('频率/hz'); ylabel('电压/v'); %-------------------------------功率谱------------------------------------ figure(6); subplot(2,1,1); youtpsd=q_d.*conj(q_d); plot(f,youtpsd); title('输出信号功率谱'); xlabel('频率/Hz'); ylabel('W/Hz'); subplot(2,1,2); ypsd=h_d.*conj(h_d); plot(f,ypsd); title('输入信号功率谱'); xlabel('频率/Hz'); ylabel('W/Hz'); %-------------------------------自相关函数-------------------------------- figure(7) subplot(2,1,1); [Rx,maxlags]=xcorr(y,'unbiased'); %信号的自相关 if fs>10000 %调整时间轴单位及标签,便于观测波形 plot(maxlags/fs*1000,Rx/max(Rx)); else plot(maxlags/fs,Rx/max(Rx)); end title('输入信号自相关函数'); xlabel('时间/s'); ylabel('R(t)'); subplot(2,1,2); [Rx1,maxlags1]=xcorr(yout,'unbiased'); %信号的自相关 if fs>10000 %调整时间轴单位及标签,便于观测波形 plot(maxlags1/fs*1000,Rx1/max(Rx)); else plot(maxlags1/fs,Rx1/max(Rx)); end title('输出信号自相关函数'); xlabel('时间/s'); ylabel('R(t)');高通滤波器如下所示：
%功能：一阶RC滤波器仿真 %说明： %1、分析了一阶RC滤波器的幅值衰减特性和相移特性 %2、分析了一阶RC滤波器的频域特性 %3、使用lsim对系统进行仿真 %4、使用FFT对原始输入信号和滤波器输出信号进行分析 %传递函数：sys=s/(s+1/RC) %========================================================================== close all; clear all; %% %输入信号的产生 %输入信号产生 fs=5000000; % fs 采样率 F=10000; % F 频率,小于采样率的一半（奈奎斯特） A=1; % A 幅度值 N=5000; % N 采样个数 %采样频率与采样个数没有什么必然关系，采样个数乘上采样频率的倒数（时间间隔）即为总的时间长度 dt=1/fs; %时间间隔 t=0:dt:(N-1)*dt; %时间向量 freqPixel=fs/N; %频率分辨率，即点与点之间频率单位 dataSourceType=2; %波形类型 0表示正弦波，1表示双音正弦波，2表示方波 F1 = 20e+3; F2 = 5e+3; F3 = 25e+3; switch dataSourceType case 0 y=A*sin(2*pi*F1*t); %频率为20khz正弦波 case 1 y=A*(sin(2*pi*F2*t)+sin(2*pi*F3*t) ); %频率为25khz和5khz双音正弦波 case 2 y=A*square(2*pi*F2*t,50);%频率为5KHz的方波 otherwise end % Tlabel='时间/s'; % if fs>10000 %调整时间轴单位及标签,便于观测波形 % t=t*1000; % Tlabel='时间/ms'; % end figure(1); subplot(2,1,1); plot(t,y); title('信号的时域波形');%对原始图像进行时域画图 xlabel('时间/s'); ylabel('电压/v'); h=fft(y,5000);%快速傅里叶变换 h_d=abs(fftshift(h));%使频域图像中间为零 f=(-N/2:N/2-1)*fs/N;%将取得时间上的点转化为频率上的点 subplot(2,1,2); plot(f,h_d/5000);%画原始图像频域上图 title('信号的频域波形'); xlabel('频率/hz'); ylabel('电压/v'); % axis([-5e+5 5e+5 0 0.6]); axis([-5e+5 5e+5 0 0.8]); %% 滤波器的设置 % f=1:1:F;%频率序列 w=2*pi*f; R=1.0010e+05;%电阻值 C=100e-12;%电容值 Fc=1/(2*pi*R*C);%截止频率 Para=R*C*1i; for n=1:length(f) A(n)=abs(1/(1+(1/(Para*w(n)))));%幅值衰减系数 P(n)=angle(1/(1+(1/(Para*w(n))))) * 180 / pi;%相移系数 end figure(2); subplot(2,1,1); plot(f,A,'r');%幅值曲线 axis([0 3e+5 0 1]); title('幅值衰减特性'); subplot(2,1,2); plot(f,P,'blue');%相位曲线 axis([0 3e+5 0 100]); title('相位特性'); Func=tf([1 0],[1 1/(R*C)]);%系统的传递函数 figure(3); bode(Func);%系统的波特图 title('幅频特性'); %% 输出信号 [yout,tout] = lsim(Func,y,t);%滤波后信号图像 %-----------------------时域图像------------------------------------ figure(4); subplot(2,1,1); plot(t,y); title('原始信号'); xlabel('时间/s'); ylabel('电压/v'); subplot(2,1,2); plot(tout,yout); title('滤波后的时域波形'); xlabel('时间/s'); ylabel('电压/v'); %------------------------------频域图像--------------------------------- q=fft(yout); q_d=abs(fftshift(q)); figure(5); subplot(2,1,1); plot(f,h_d/5000); title('输入信号的频谱'); xlabel('频率/hz'); ylabel('电压/v'); % axis([-5e+5 5e+5 0 0.6]); axis([-5e+5 5e+5 0 0.8]); subplot(2,1,2); plot(f,q_d/5000); title('滤波后的频谱'); xlabel('频率/hz'); ylabel('电压/v'); % axis([-5e+5 5e+5 0 0.6]); axis([-5e+5 5e+5 0 0.8]); %-------------------------------功率谱------------------------------------ figure(6); subplot(2,1,1); ypsd=h_d.*conj(h_d); plot(f,ypsd); title('输入信号功率谱'); xlabel('频率/Hz'); ylabel('W/Hz'); % axis([-5e+5 5e+5 0 0.6]); axis([-5e+5 5e+5 0 12e+6]); subplot(2,1,2); youtpsd=q_d.*conj(q_d); plot(f,youtpsd); title('输出信号功率谱'); xlabel('频率/Hz'); ylabel('W/Hz'); % axis([-5e+5 5e+5 0 0.6]); axis([-5e+5 5e+5 0 5e+6]); %-------------------------------自相关函数-------------------------------- figure(7) subplot(2,1,1); [Rx,maxlags]=xcorr(y,'unbiased'); %信号的自相关 if fs>10000 %调整时间轴单位及标签,便于观测波形 plot(maxlags/fs*1000,Rx/max(Rx)); else plot(maxlags/fs,Rx/max(Rx)); end title('输入信号自相关函数'); xlabel('时间/s'); ylabel('R(t)'); subplot(2,1,2); [Rx1,maxlags1]=xcorr(yout,'unbiased'); %信号的自相关 if fs>10000 %调整时间轴单位及标签,便于观测波形 plot(maxlags1/fs*1000,Rx1/max(Rx)); else plot(maxlags1/fs,Rx1/max(Rx)); end title('输出信号自相关函数'); xlabel('时间/s'); ylabel('R(t)');注意传递函数tf的用法。
tf([1 1],[1 1 1])表示响应函数(s+1)/(s^2 + s + 1)
