【公開課】數(shù)字信號分析理論與實(shí)踐（基于matlab） - 華中科技大學(xué)（數(shù)字信
課程中的matlib版本較低，很多函數(shù)發(fā)生了變化。在這里補(bǔ)充一下課程中的實(shí)例程序。
使用到的控件：按鈕、編輯字段（數(shù)值）、滑塊、坐標(biāo)區(qū)。
需要注意的是：最后老師演示的實(shí)例，滑動滑塊后自動顯示了波形，這是因?yàn)槔蠋煂⒊绦蛱砑釉诹嘶瑝K的回調(diào)函數(shù)中。所使用的程序不需要額外改變。
%在標(biāo)尺中可以設(shè)置坐標(biāo)軸的范圍來實(shí)現(xiàn)對整體信號的截取，程序中整個聲音時長為1s，通過坐標(biāo)軸范圍來顯示0.01s中的波形
% 生成音頻文件需要為全局變量，通過使用global關(guān)鍵詞來設(shè)置
??????global player;
??????F=app.Slider.Value;
??????Fs=44100;
??????dt=1.0/Fs;
??????T=1;
??????N=T/dt;
??????t=linspace(0,T,N);
??????x=0.3*sin(2*pi*F*t);
??????plot(app.UIAxes,t,x);
??????app.EditField.Value = F;
??????player = audioplayer(x,Fs);%生成音頻文件
??????play(player);%播放音頻文件
P16 作業(yè)-----------------------------------------------
%用11025Hz的采樣頻率，用Matlab生成時間長度為1s的500Hz的正弦波
%方波等標(biāo)準(zhǔn)信號，并用plot函數(shù)繪制信號波形，用sound函數(shù)播放信號聲音
Fs=11025;
dt=1.0/Fs;
T=1;
N=T/dt;
t=(0:N-1)/N;
%正弦波
%x=0.3*sin(2*pi*500*t);
%方波
%x=0.3*square(2*pi*500*t);
%鋸齒波
%x=0.3*sawtooth(2*pi*500*t,1/2);
%白噪聲
%x=0.3*randn(1,N);
plot(t,x);
axis([0,0.01,-0.5,0.5]);
player = audioplayer(x,Fs);%生成音頻文件
play(player);%播放音頻文件
P19 電子琴加泛音-----------------------------------
只以1號鍵為例
global player
??????A=linspace(0,0.9,8800);
??????D=linspace(0.9,0.8,2200);
??????S=linspace(0.8,0.8,18000);
??????R=linspace(0.8,0,15100);
??????global adsr
??????adsr=[A,D,S,R];
??????Fs= 44100;
??????dt=1.0/Fs;
??????T=1;
??????N=T/dt;
??????t=linspace(0,T,N);
??????x=0.3*sin(2*pi*247*t);
??????y=x.*adsr;
??????plot(app.UIAxes,t,y);
??????player = audioplayer(y,Fs);
??????play(player);
項(xiàng)目1：數(shù)字信號發(fā)生器
效果圖：?
回調(diào)函數(shù)列表：
?控件列表：
代碼如下：
??properties (Access = private)%私有屬性
????Timer_id;
????type;
??end
???
??methods (Access = private)%私有函數(shù)
?????
????function my_callback_fcn(app)
??????global Fs
??????global x
??????global N
??????global F
??????global A
??????global z
??????Fs = 44100;
??????dt=1.0/Fs;
??????T=1;
??????N=T/dt;
??????t=linspace(0,T,N);
??????F=app.FreSlider.Value;
??????A=app.AmpSlider.Value;
??????if app.type==1
????????x=0.5*A*randn(1,N);
??????end
??????if app.type==2
????????x=A*sin(2*pi*F*t+z);
??????end
??????if app.type==3
????????x=A*square(2*pi*F*t+z);
??????end
??????if app.type==4
????????x=A*sawtooth(2*pi*F*t+z);
??????end
??????z=z+0.25;
??????plot(app.UIAxes,t,x);
??????grid(app.UIAxes);
???????
????end
????%定時器初始化，重復(fù)模式
????function timer_init(app)
??????app.Timer_id=timer;
??????app.Timer_id.Period=1.0; %周期
??????app.Timer_id.ExecutionMode='fixedRate';
??????app.Timer_id.TasksToExecute=Inf;
??????app.Timer_id.TimerFcn=@(~,~) my_callback_fcn(app);
????end
????%定時器啟動
????function timer_start(app)
??????start(app.Timer_id);
????end
????%定時停止
????function timer_stop(app)
??????stop(app.Timer_id);
????end
????%刪除定時器
????function timer_delete(app)
??????delete(app.Timer_id);
????end
??end
???
??% Callbacks that handle component events
??methods (Access = private)
????% Value changed function: FreSlider
????function FreSliderValueChanged(app, event)
??????app.EditFreEditField.Value=event.Value;??????
????end
????% Value changed function: AmpSlider
????function AmpSliderValueChanged(app, event)
??????app.EditAmpEditField.Value=event.Value;???????
????end
????% Value changed function: EditFreEditField
????function EditFreEditFieldValueChanged(app, event)
??????app.FreSlider.Value=event.Value;??????
????end
????% Value changed function: EditAmpEditField
????function EditAmpEditFieldValueChanged(app, event)
??????app.AmpSlider.Value=event.Value;?????
????end
????% Button pushed function: PlayButton
????function PlayButtonPushed(app, event)
??????global player
??????global Fs
??????global x
??????y=x*0.01;
??????player=audioplayer(y,Fs);
??????play(player);
????end
????% Button pushed function: NoiseButton
????function NoiseButtonPushed(app, event)
??????app.type=1;
??????my_callback_fcn(app);
????end
????% Button pushed function: SineButton
????function SineButtonPushed(app, event)
??????app.type=2;
??????my_callback_fcn(app);
????end
????% Button pushed function: SquareButton
????function SquareButtonPushed(app, event)
??????app.type=3;
??????my_callback_fcn(app);
????end
????% Button pushed function: TrangleButton
????function TrangleButtonPushed(app, event)
??????app.type=4;
??????my_callback_fcn(app);
????end
????% Button pushed function: RunButton
????function RunButtonPushed(app, event)
??????timer_init(app);
??????timer_start(app);
????end
????% Button pushed function: StopButton
????function StopButtonPushed(app, event)
??????ts=timerfind;
??????if ~isempty(ts)
????????timer_stop(app);
????????timer_delete(app);
??????end?????
????end
??end
P24- 作業(yè)----------------------------------------------
本人將兩個作業(yè)合并為一個。
效果圖：
所使用到的控件：
程序代碼如下：
???
??properties (Access = private)
????x % Description
????dt
????N
????t
????Fs=5120;
????type_num=1;
??end
???
??methods (Access = private)
?????
????function wave_sin(app)
??????pinlv=app.Slider.Value;
??????fuzhi=app.Slider_2.Value;
??????app.dt=1.0/app.Fs;
??????T=1;
??????app.N=T/app.dt;
??????app.t=(0:app.N-1)/app.N*T;
??????app.x=fuzhi*sin(2*pi*pinlv*app.t);
??????jisuan(app);
????end??
?????
????function daoshu(app)%對信號先微分再積分
??????y1=app.x;
??????y1(1)=app.x(1);
??????for k=2:1:app.N-1
????????y1(k)=(app.x(k+1) - app.x(k-1))/(2*app.dt);
??????end
??????y1(app.N)=app.x(app.N);
??????plot(app.UIAxes2,app.t,y1);
??????y2=y1;
??????y2(1)=0;
??????for k=2:1:app.N-1
????????y2(k)=y2(k-1) + (app.dt*(y1(k)+y1(k-1)))/2;
??????end
??????plot(app.UIAxes3,app.t,y2);
????end
?????
????function jisuan(app)
???????%分析幅值
??????fz=max(app.x);
??????app.EditField.Value=fz;
??????%分析峰峰值
??????ffz=max(app.x)-min(app.x);
??????app.EditField_2.Value=ffz;
??????%分析平均值
??????pjz=mean(app.x);
??????app.EditField_3.Value=pjz;
??????%分析有效值
??????cz=app.x-pjz;%波形數(shù)組每一個元素與平均值的差值
??????czpf=cz.^2;%差值平方
??????czpfpj=sum(czpf)/app.N;%差值平方求和取平均
??????yxz=sqrt(czpfpj);%開方求得有效值
??????app.EditField_5.Value=yxz;
??????%標(biāo)準(zhǔn)差
??????app.EditField_6.Value=std(app.x);
??????%頻率和初始相位
??????p=max(app.x);q=min(app.x);n=1;ti=(0);
??????at=0.8*(p-q)+q;
??????for k=2:1:app.N-1
????????if(app.x(k-1)<at && app.x(k)<=at && app.x(k+1)>at && app.x(k+2)>at)
??????????ti(n) = k;
??????????n = n+1;
????????end
??????end
??????T=((ti(2) - ti(1)))*app.dt;
??????F=1.0/T;
??????Q=(T-ti(1)*app.dt);%計(jì)算出初始相位
??????app.EditField_9.Value=Q;
??????app.EditField_4.Value=F;
??????plot(app.UIAxes,app.t,app.x);??????
????end
?????
????function wave_square(app)
??????pinlv=app.Slider.Value;
??????fuzhi=app.Slider_2.Value;
??????app.dt=1.0/app.Fs;
??????T=1;
??????app.N=T/app.dt;
??????app.t=(0:app.N-1)/app.N*T;
??????app.x=fuzhi*square(2*pi*pinlv*app.t);
??????jisuan(app);
????end
?????
????function wave_sawtooth(app)
??????pinlv=app.Slider.Value;
??????fuzhi=app.Slider_2.Value;
??????app.dt=1.0/app.Fs;
??????T=1;
??????app.N=T/app.dt;
??????app.t=(0:app.N-1)/app.N*T;
??????app.x=fuzhi*sawtooth(2*pi*pinlv*app.t);
??????jisuan(app);????
????end
?????
????function select(app)
??????switch (app.type_num)
????????case 1
??????????try??????
????????????wave_sin(app);
????????????daoshu(app);
??????????catch ErrorInfo
????????????warndlg('數(shù)據(jù)輸入有誤?。?#39;);
??????????end
????????case 2
??????????try????????????
????????????wave_square(app);
????????????daoshu(app);
??????????catch ErrorInfo
????????????warndlg('數(shù)據(jù)輸入有誤??！');
??????????end
????????case 3
??????????try
????????????wave_sawtooth(app);
????????????daoshu(app);
??????????catch ErrorInfo
????????????warndlg('數(shù)據(jù)輸入有誤??！');
??????????end
??????end
????end
??end
??% Callbacks that handle component events
??methods (Access = private)
????% Value changed function: Slider_2
????function Slider_2ValueChanged(app, event)
??????app.EditField_7.Value = app.Slider_2.Value;
??????select(app);
????end
????% Value changed function: EditField_7
????function EditField_7ValueChanged(app, event)
??????app.Slider_2.Value = app.EditField_7.Value;
??????select(app);
????end
????% Value changed function: Slider
????function SliderValueChanged(app, event)
??????app.EditField_8.Value = app.Slider.Value;
??????select(app);
????end
????% Value changed function: EditField_8
????function EditField_8ValueChanged(app, event)
??????app.Slider.Value = app.EditField_8.Value;
??????select(app);
????end
????% Button pushed function: Button
????function ButtonPushed(app, event)
??????try?
????????app.type_num=1;
????????wave_sin(app);
????????daoshu(app);
??????catch ErrorInfo
????????warndlg('數(shù)據(jù)輸入有誤?。?#39;);
??????end
????end
????% Button pushed function: Button_2
????function Button_2Pushed(app, event)
??????try?
????????app.type_num=2;
????????wave_square(app);
????????daoshu(app);
??????catch ErrorInfo
????????warndlg('數(shù)據(jù)輸入有誤??！');
??????end
????end
????% Button pushed function: Button_3
????function Button_3Pushed(app, event)
??????try?
????????app.type_num=3;
????????wave_sawtooth(app);
????????daoshu(app);
??????catch ErrorInfo
????????warndlg('數(shù)據(jù)輸入有誤?。?#39;);
??????end
????end
??end
P27---------------------------------------------------
Fs=5120;N=1024;
dt=1.0/Fs;T=dt*N;
t=linspace(0,T,N);
x=10*sin(2*pi*100*t)+10/3*sin(2*pi*3*100*t);
plot(t,x);
title('原始信號');
y=fft(x,N);
a=real(y);b=imag(y);
figure;%創(chuàng)建圖窗窗口
subplot(2,1,1);
plot(a);
title('實(shí)頻譜');
subplot(2,1,2);
plot(b);
title('虛頻譜');
f=linspace(0,Fs/2,N/2);%x坐標(biāo)換為頻率
A1=abs(y)/(N/2);%幅值量綱還原
Q1=angle(y)*180/pi;
figure;%創(chuàng)建圖窗窗口
subplot(2,1,1);
plot(f,A1(1:N/2));%不顯示負(fù)頻率部分
title('幅頻譜');
subplot(2,1,2);
plot(f,Q1(1:N/2));%不顯示負(fù)頻率部分
title('相頻譜');
A2=A1.^2;%功率譜
P2=20*log10(A2);%對數(shù)功率譜
figure;
subplot(2,1,1);
plot(f,A2(1:N/2));
title('功率譜');
subplot(2,1,2);
plot(f,P2(1:N/2));
title('對數(shù)功率譜');
P27 作業(yè)---------------------------------------------
Fs=11025;N=1024;
dt=1.0/Fs;T=dt*N;
t=linspace(0,T,N);
%正弦波
x=10*sin(2*pi*200*t);
y=fft(x,N);
f=linspace(0,Fs/2,N/2);%x坐標(biāo)換為頻率
A1=abs(y)/(N/2);%幅值量綱還原
subplot(2,1,1);
plot(t,x);
title('正弦波');
subplot(2,1,2);
plot(f,A1(1:N/2));%不顯示負(fù)頻率部分
title('幅頻譜');
%方波
x1=10*square(2*pi*200*t);
y1=fft(x1,N);
A2=abs(y1)/(N/2);%幅值量綱還原
figure;
subplot(2,1,1);
plot(t,x1);
title('方波');
subplot(2,1,2);
plot(f,A2(1:N/2));%不顯示負(fù)頻率部分
title('幅頻譜');
%三角波
x2=10*sawtooth(2*pi*200*t);
y2=fft(x2,N);
A3=abs(y2)/(N/2);%幅值量綱還原
figure;
subplot(2,1,1);
plot(t,x2);
title('三角波');
subplot(2,1,2);
plot(f,A3(1:N/2));%不顯示負(fù)頻率部分
title('幅頻譜');
P34------------------------------------------------------
【備注】：只是完成了課程中老師所提及到的功能，更多細(xì)節(jié)可以再優(yōu)化優(yōu)化，就讓各位讀者自行探索。
效果圖：
程序：
???
??properties (Access = private)
????Timer_id
????wave
????num=0;
????N=24810;?% 24810 = length(doubleArray)/2
????Z1=0.1;
????Z2=0.1;
??end
???
??methods (Access = private)
????%音頻處理
????function Recording(app)
??????recObj = audiorecorder(11025,16,1);
??????recordblocking(recObj,5);
??????doubleArray = getaudiodata(recObj);
??????doubleArray = doubleArray(5506:end);%去除錄音數(shù)據(jù)頭
??????ylim(app.UIAxes,[-app.Z1 app.Z1]);
??????plot(app.UIAxes,doubleArray);
??????switch app.num
????????case 0
??????????pinlv(app,doubleArray);
????????case 1
??????????%注意兩個數(shù)組的大小
??????????H=2*hamming(app.N*2).*doubleArray;%2為修正系數(shù)
??????????pinlv(app,H);
????????case 2
??????????B=2.381*blackman(app.N*2).*doubleArray;
??????????pinlv(app,B);
????????case 3
??????????P=4.545*flattopwin(app.N*2).*doubleArray;
??????????pinlv(app,P);
??????end
????end
?????
????%fft
????function pinlv(app,wave)??????
??????y_fft=fft(wave);
??????f=linspace(0,5505,app.N);
??????A=abs(y_fft)/app.N;
??????ylim(app.UIAxes2,[0 0.01*app.Z2]);
??????plot(app.UIAxes2,f,A(1:app.N));??????
????end
?????
????%定時器初始化，重復(fù)模式
????function timer_init(app)
??????app.Timer_id=timer;
??????app.Timer_id.Period=5.2; %周期
??????app.Timer_id.ExecutionMode='fixedRate';
??????app.Timer_id.TasksToExecute=Inf;
??????app.Timer_id.TimerFcn=@(~,~) Recording(app);
????end
?????
????%定時器啟動
????function timer_start(app)
??????start(app.Timer_id);
????end
?????
????%定時停止
????function timer_stop(app)
??????stop(app.Timer_id);
????end
?????
????%刪除定時器
????function timer_delete(app)
??????delete(app.Timer_id);
????end
??end
???
??% Callbacks that handle component events
??methods (Access = private)
????% Button pushed function: Button
????function ButtonPushed(app, event)
??????timer_init(app);
??????timer_start(app);
????end
????% Button pushed function: Button_2
????function Button_2Pushed(app, event)
??????timer_delete(app);
????end
????% Button pushed function: Button_3
????function Button_3Pushed(app, event)
??????app.num = 1;
????end
????% Button pushed function: Button_4
????function Button_4Pushed(app, event)
??????app.num = 2;
????end
????% Button pushed function: Button_5
????function Button_5Pushed(app, event)
??????app.num = 3;
????end
????% Value changed function: Slider
????function SliderValueChanged(app, event)
??????app.Z1 = 0.1*2*app.Slider.Value;???
????end
????% Value changed function: Slider_2
????function Slider_2ValueChanged(app, event)
??????app.Z2 = 0.1*2*app.Slider_2.Value;
????end
??end
P36 第五章-----------------------------------------
信號的時差域相關(guān)分析
信號的相關(guān)分析是一種分析兩個信號之間或一個信號自身的時間依存關(guān)系和相似程度的方法。
工程上，人們關(guān)心的是信號不同時刻的相似程度，但不太關(guān)心其具體值，這時相關(guān)函數(shù)可簡化為：
x=y 自相關(guān)；x!=y 互相關(guān)
數(shù)字信號離散計(jì)算公式：
聲音素材網(wǎng)址：https://sc.chinaz.com/yinxiao/220314192200.htm
程序：
[y,Fs] = audioread('beat.wav');%讀取聲音文件，需要在同一文件夾下
y1=y(:,1);%取單聲道
subplot(2,1,1);
plot(y1);
s=xcorr(y1,'unbiased');%自相關(guān)函數(shù)
% ‘unbiased' 避免重疊失真
subplot(2,1,2);
plot(s);
結(jié)果：下圖中相關(guān)函數(shù)圖形，橫坐標(biāo)中點(diǎn)為零點(diǎn)。
相關(guān)函數(shù)的性質(zhì)：
（1）自相關(guān)函數(shù)是偶函數(shù)
（2）當(dāng)t=0時，自相關(guān)函數(shù)具有最大值
（3）周期信號的自相關(guān)函數(shù)仍然是同頻率的周期信號，但不保留原信號的相位信息
（4）隨機(jī)噪聲信號的自相關(guān)函數(shù)將隨t的增大快速衰減
（5）兩周期信號的互相關(guān)函數(shù)仍然是同頻率的周期信號，且保留原有信號的相位信息
（6）兩個非同頻率的周期信號互不相關(guān)
P40----------------------------------------------------
在數(shù)學(xué)中，概率密度函數(shù)(Probability Density Function)是一個描述信號的取值在某個確定的取值點(diǎn)附近的可能性的函數(shù)，定義為概率分布函數(shù)的倒數(shù)。
以幅值大小為橫坐標(biāo)，以每個幅值間隔內(nèi)出現(xiàn)的概率為縱坐標(biāo)進(jìn)行統(tǒng)計(jì)分析的方法。它反映了信號落在不同幅值強(qiáng)度區(qū)域內(nèi)的概率情況。
概率分布函數(shù)
概率分布函數(shù)(Cumulative distribution function)是信號幅值小于或等于某值R的概率，其定義為：
概率分布函數(shù)又稱之為累積概率，表示了落在某一區(qū)間的概率。
直方圖
以幅值為橫坐標(biāo)，以每個幅值間隔內(nèi)出現(xiàn)的頻次為縱坐標(biāo)進(jìn)行統(tǒng)計(jì)分析，稱為直方圖。
直方圖歸一化后可轉(zhuǎn)換為概率密度函數(shù)。
直方圖均衡，主要作用是提高圖像的清晰度和紋理的識別度，可能會使得色彩失真，更多的是用在灰度圖像的處理。
P42------------------------------------------------------
效果圖：
代碼：
function ButtonValueChanged(app, event)
??????[filename, filepath]=uigetfile({'*.bmp;*.jpg;*.png'},'選擇圖片文件'); %選擇文件
??????if isequal(filename,0) ||?isequal(filepath,0)%判斷是否選中圖片
????????errordlg('沒有選中文件，請重新選擇','出錯');
????????return;
??????else
?????????set(app.EditField,'Value',[filepath, filename]);%獲取文件的路徑和名字
??????end
??????road=get(app.EditField,'Value');
??????set(gcbf,'UserData',road);% GUI中參數(shù)的傳遞，詳情https://www.cnblogs.com/jmliao/p/5628521.html
????end
????% Value changed function: Button_2
????function Button_2ValueChanged(app, event)
?????road=get(gcbf,'UserData');
?????app.Image.ImageSource=road;%顯示原圖像
?????O_image=imread(road);%讀入原始圖像
?????Image_gray=im2gray(O_image);%轉(zhuǎn)換為灰度圖像
?????imwrite(Image_gray,'myGray.png');%將灰度圖像保存，圖像控件只能顯示真色彩圖像，暫未找到直接顯示灰度圖像的方法
?????app.Image_3.ImageSource='myGray.png';
?????Image_A_gray=histeq(Image_gray);%直方圖均衡強(qiáng)化對比度
?????imwrite(Image_A_gray,'myHisteq.png');
?????app.Image_4.ImageSource='myHisteq.png';
????end
P43------------------------------------------------
白噪聲信號的幅值域分析
效果圖：
程序：
x=-4:0.02:4;
y=randn(10000,1);
subplot(4,1,1);
plot(y);
[N, edges] = histcounts(y,x);%N-bin中的計(jì)數(shù)，edges-bin的邊界
centers = (edges(1:end-1) + edges(2:end))./2;%取邊界中心值
subplot(4,1,2);stem(centers,N); %繪制離散序列數(shù)據(jù)
pdf=N/length(y);%歸一化,也可以直接使用示例程序中的方式進(jìn)行歸一化
subplot(4,1,3);plot(centers,pdf);
subplot(4,1,4);cdfplot(y);%累計(jì)分布函數(shù)圖像
補(bǔ)充：(可以參考一下，便于理解直方圖的參數(shù)和屬性）
示例程序（來源網(wǎng)絡(luò)）：
x = randn(1,1000); %Generate some data to plot
figure(1); clf;
subplot(2,1,1); hold on;
h = histogram(x,'normalization','probability');
title('Plotted as a histogram');
ylabel('Frequency');
subplot(2,1,2); hold on;
[N, edges] = histcounts(x,'normalization','probability');
centers = (edges(1:end-1) + edges(2:end))./2; %histcounts gives the bin edges, but we want to plot the bin centers
plot(centers,N,'ko');
title('Plotted as points');
ylabel('Frequency');
自己畫的，便于理解：
直方圖，概率質(zhì)量函數(shù)和概率密度函數(shù)
https://blog.csdn.net/lovehua365/article/details/79058815
頻譜、能譜、功率譜、倍頻程譜、1/3 倍頻程譜
https://blog.csdn.net/liyuanbhu/article/details/42675765
2倍頻程、1倍頻程、1/3倍頻程介紹
https://blog.csdn.net/chenhuanqiangnihao/article/details/122079324
P46 信號的數(shù)字濾波-------------------------------
測量中除傳感器信號外，干擾也會出現(xiàn)，他們與信號疊加在一起，扭曲測量結(jié)果。濾波器是一種選頻裝置，可以使信號中特定頻率成分通過，而衰減其他頻率成分。
濾波有一個前提條件，信號的頻率是不重疊的。
P47---------------------------------------------
數(shù)字濾波和模擬濾波的區(qū)別
濾波位置的不同；采樣頻率不同；有效A/D位數(shù)；
當(dāng)干擾信號頻率較高時，應(yīng)優(yōu)先模擬濾波將高頻的干擾信號濾除，這樣可以降低采樣頻率。
當(dāng)干擾信號的量程較大時，直接進(jìn)行數(shù)字濾波會使得A/D滿量程精度較小。
濾波器評價(jià)指標(biāo)：過渡區(qū)陡峭；脈沖響應(yīng)緊支集。
P48---------------------------------------------
Fs=2048;dt=1.0/Fs;T=1;
N=T/dt;t=[0:N-1]/N;
x1=sin(2*pi*50*t)+sin(2*pi*300*t)+sin(2*pi*500*t);
subplot(4,1,1);
plot(t,x1);axis([0,0.1,-2,2]);
P=fft(x1,N);
Pyy=2*sqrt(P.*conj(P))/N;%幅值量綱還原
f=linspace(0,Fs/2,N/2);
subplot(4,1,2);plot(f,Pyy(1:N/2));
for k=1:N
??P1(k)=real(P(k))+i*(imag(P(k)));
end
for k=400:N-400
??P1(k)=0;
end
for k=1:200
??P1(k)=0;
end
for k=N-200:N
??P1(k)=0;
end
Pyy=2*sqrt(P1.*conj(P1))/N;
f=linspace(0,Fs/2,N/2);
subplot(4,1,3);plot(f,Pyy(1:N/2));
x2=ifft(P1);
subplot(4,1,4);plot(t,x2);
axis([0,0.1,-2,2]);
P49---------------------------------------------
%數(shù)字差分--簡單的高通濾波器
Fs=5000;a=5;
f=2;T=1;
dt=1.0/Fs;
N=T/dt;
t=linspace(0,T,N);
y=a*sin(2*3.14*f*t)+0.3*sin(100*3.14*f*t);
plot(t,y);
x(1)=0;x(N)=0;
for i=2:N-1
??x(i)=(y(i+1)-y(i-1))/2;
end
figure;
plot(t,x);
%數(shù)據(jù)平滑--簡單的低通濾波器
Fs=5000;a=5;
f=2;T=1;
dt=1.0/Fs;
N=T/dt;
t=linspace(0,T,N);
y=a*sin(2*pi*f*t)+0.8*sin(1000*pi*f*t);
subplot(2,1,1);
plot(t,y);
x=y;
for i=1:N-5
??x(i)=[y(i)+y(i+1)+y(i+2)+y(i+3)+y(i+4)]/5;
end
subplot(2,1,2);
plot(t,x);
作業(yè)
Fs=5120;
T=1;
dt=1.0/Fs;
N=T/dt;
t=linspace(0,T,N);
y=sin(2*pi*10*t)+sin(2*pi*500*t);
subplot(5,1,1);
plot(t,y);
%5點(diǎn)平滑濾波
x=y;
for a=1:N-5
??x(a)=[y(a)+y(a+1)+y(a+2)+y(a+3)+y(a+4)]/5;
end
subplot(5,1,2);
plot(t,x);
%頻域?yàn)V波
P=fft(y,N);
Pyy=sqrt(P.*conj(P))/N;%幅值量綱還原
d=linspace(0,Fs,N);
subplot(5,1,3);
plot(d,Pyy);
for k=1:N
??P1(k)=real(P(k))+i*(imag(P(k)));
end
for k=200:N-200
??P1(k)=0;
end
Pyy1=sqrt(P1.*conj(P1))/N;
d=linspace(0,Fs,N);
subplot(5,1,4);
plot(d,Pyy1);
x2=ifft(P1);
subplot(5,1,5);plot(t,x2);
P54--------------------------------------------------------------------------
Fs=2048;dt=1/Fs;
T=1;N=T/dt;t=[0:N-1]/N;
x1=sin(2*pi*50*t)+sin(2*pi*300*t)+sin(2*pi*500*t);
subplot(3,1,1);plot(t,x1);
axis([0,0.1,-2,2]);P=fft(x1,N);
Pyy = 2*sqrt(P.*conj(P))/N;
f=linspace(0,Fs/2,N/2);
subplot(3,1,2);plot(f,Pyy(1:N/2));
b=fir1(48,0.1);
x2=filter(b,1,x1);
subplot(3,1,3);plot(t,x2);axis([0,0.1,-2,2])
figure;
c=fir1(48,[0.2,0.4]);
x3=filter(c,1,x1);
subplot(2,1,1);plot(t,x3);axis([0,0.1,-2,2])
d=fir1(48,0.4,'high');
x4=filter(d,1,x1);
subplot(2,1,2);plot(t,x4);axis([0,0.1,-2,2])
%%?原始信號
Fs=5120;dt=1/Fs;
T=1;N=T/dt;t=[0:N-1]/N;
%含有低(100Hz)、中(600Hz)、高(1000Hz)三個不同頻率成分的測試信號
x1=sin(2*pi*100*t)+sin(2*pi*600*t)+sin(2*pi*1000*t);
subplot(2,1,1);plot(t,x1);
axis([0,0.1,-3,3]);P=fft(x1,N);
Pyy = 2*sqrt(P.*conj(P))/N;
f=linspace(0,Fs/2,N/2);
subplot(2,1,2);plot(f,Pyy(1:N/2));
%% fir1
figure;
b=fir1(48,0.1);
%48-代表濾波器的階數(shù)，指在濾波器的傳遞函數(shù)中有幾個極點(diǎn)，同時決定了轉(zhuǎn)折區(qū)的下降速度
x2=filter(b,1,x1);
subplot(3,1,1);plot(t,x2);axis([0,0.1,-2,2])
c=fir1(48,[0.1,0.3]);
x3=filter(c,1,x1);
subplot(3,1,2);plot(t,x3);axis([0,0.1,-2,2])
d=fir1(48,0.3,'high');
x4=filter(d,1,x1);
subplot(3,1,3);plot(t,x4);axis([0,0.1,-2,2])
%% fir2
figure
f1=[0 0.1 0.1 1];%低通
mhi1=[1 1 0 0];
b1=fir2(48,f1,mhi1);
x2_fir2=filter(b1,1,x1);
subplot(3,1,1);plot(t,x2_fir2);axis([0,0.1,-2,2])
f2=[0 0.1 0.1 0.3 0.3 1];
mhi2=[0 0 1 1 0 0];
c1=fir2(48,f2,mhi2);
x3_fir2=filter(c1,1,x1);
subplot(3,1,2);plot(t,x3_fir2);axis([0,0.1,-2,2])
f3=[0 0.3 0.3 1];
mhi3=[0 0 1 1];
d1=fir2(48,f3,mhi3);
x4_fir2=filter(d1,1,x1);
subplot(3,1,3);plot(t,x4_fir2);axis([0,0.1,-2,2])
%% firls
figure
% f1=[0 0.1 0.1 1];%低通
% mhi1=[1 1 0 0];
b2=firls(48,f1,mhi1);
x2_firls=filter(b2,1,x1);
subplot(3,1,1);plot(t,x2_firls);axis([0,0.1,-2,2])
% f2=[0 0.1 0.1 0.3 0.3 1];
% mhi2=[0 0 1 1 0 0];
c2=firls(48,f2,mhi2);
x3_firls=filter(c2,1,x1);
subplot(3,1,2);plot(t,x3_firls);axis([0,0.1,-2,2])
% f3=[0 0.3 0.3 1];
% mhi3=[0 0 1 1];
d2=firls(48,f3,mhi3);
x4_firls=filter(d2,1,x1);
subplot(3,1,3);plot(t,x4_firls);axis([0,0.1,-2,2])
