原始信號由20Hz、50Hz和100Hz三種頻率的正弦波組成，如上圖所示。
在Figure圖上畫框截取待分析的數據，如下所示：
所得結果如下所示，
從上圖中可以清晰的看出有三種頻率，而且從時頻圖中可以看到三種頻率發生在不 同時刻。
點擊單邊傅里葉頻譜，放大頻譜如下：
畫框可以選定濾波范圍，如帶通濾波，選擇如下框，可得結果：
MATLAB中具體代碼如下：
function frequency_analysis
clc close all
Ts = 0.001;
Fs = 1/Ts;
f1 = 20;
f2 = 50;
f3 = 100;
dt = 0.2;
t1 = (0:Ts:dt-Ts) + 0; t2 = (0:Ts:dt-Ts) + dt;
t3 = (0:Ts:dt-Ts) + 2*dt;
y1 = sin(2*pi*f1*t1);
y2 = sin(2*pi*f2*t2);
y3 = sin(2*pi*f3*t3);
t = [t1 t2 t3];
y = [y1 y2 y3];
figure
plot(t,y)
xlim([t(1) t(end)])
ylim([min(y) max(y)])
xlabel('時間t')
ylabel('信號y(t)')
title('原始信號')
%
% 以下為標準化程序
% 上面的Ts和Fs要給定，后面會用到
set(gcf,'WindowButtonDownFcn',@BtndownFcn);
function BtndownFcn(h,evt)
temppt = get(gca,'CurrentPoint');
startpt.x = temppt(1,1);
startpt.y = temppt(1,2);
endpt = startpt;
height = 0.0001;
width = 0.0001;
rectangle('Position',[startpt.x, startpt.y, width, height],'Tag','rangerect'); set(h,'WindowButtonMotionFcn',@BtnmoveFcn);
set(h,'WindowButtonUpFcn',@BtnupFcn);
function BtnmoveFcn(h,evt)
temppt = get(gca,'CurrentPoint');
endpt.x = temppt(1,1);
endpt.y = temppt(1,2);
width = abs(endpt.x-startpt.x)+0.00001;
height = abs(endpt.y-startpt.y)+0.00001;
hrect = findobj('Tag','rangerect');
set(hrect,'Position',[min(startpt.x, endpt.x), min(startpt.y, endpt.y), width, height]);
end
function BtnupFcn(h,evt)
set(h,'WindowButtonMotionFcn','');
set(h,'WindowButtonUpFcn','');
hrect = findobj('Tag','rangerect');
delete(hrect);
BtnUp_Spectrum_Analysis(startpt.x, endpt.x)
end
% 頻譜分析
function BtnUp_Spectrum_Analysis(starttime, endtime)
hline = findobj(gca,'type','line');
time = get(hline,'xdata');
laser = get(hline,'ydata');
% 得到截取的分析數據
tempx = find(time >= min(starttime, endtime));
startindex = tempx(1);
tempx = find(time <= max(starttime, endtime));
endindex = tempx(end);
yt = laser(startindex:endindex);
yt = ytmean(yt);
t = time(startindex:endindex);
% 傅里葉變換
[Yf, f] = Spectrum_Calc(yt,Fs);
% 小波變換
scale = 1:50;
cw2 = cwt(yt,scale,'morl');
% 作圖
figure
subplot(231)
% 截取的頻譜分析的數據
plot(t, yt)
xlim([t(1) t(end)])
ylim([min(yt),max(yt)]);
title('頻譜分析數據')
xlabel('時間t')
ylabel('截取的數據y(t)')
h1 = subplot(234);
% 單邊傅里葉變換分析頻譜
plot(f,Yf,'-')
title('單邊傅里葉頻譜')
xlabel('頻率Hz')
ylabel('|Y(f)|')
set(h1,'ButtonDownFcn',@BtnDown_Filter_fcn)
function BtnDown_Filter_fcn(h,evt)
fig_fft = figure;
plot(f,Yf,'-')
title('單邊傅里葉頻譜')
xlabel('頻率Hz')
ylabel('|Y(f)|')
xlim([0 Fs/2])
ylim([0 max(Yf)])
% 構造右鍵菜單
filter_flag = 1;
hcmenu = uicontextmenu;
uimenu(hcmenu, 'Label', '低通濾波', 'Callback', @hcb1);
uimenu(hcmenu, 'Label', '高通濾波', 'Callback', @hcb2);
uimenu(hcmenu, 'Label', '帶通濾波', 'Callback', @hcb3);
uimenu(hcmenu, 'Label', '帶阻濾波', 'Callback', @hcb4);
set(gca,'UIContextMenu',hcmenu);
function hcb1(h,evt)
filter_flag = 1;
end
function hcb2(h,evt)
filter_flag = 2;
end
function hcb3(h,evt)
filter_flag = 3;
end
function hcb4(h,evt)
filter_flag = 4;
end
set(fig_fft,'WindowButtonDownFcn',@BtndownFcn);
function BtndownFcn(h,evt)
if
strcmp(get(h,'SelectionType'),'normal')
temppt = get(gca,'CurrentPoint');
startpt.x = temppt(1,1);
startpt.y = temppt(1,2);
endpt = startpt;
height = 0.0001;
width = 0.0001;
rectangle('Position',
[startpt.x, startpt.y, width, height],'Tag','rangerect_fft'); set(h,'WindowButtonMotionFcn',@BtnmoveFcn);
set(h,'WindowButtonUpFcn',@BtnupFcn);
end
function BtnmoveFcn(h,evt)
temppt = get(gca,'CurrentPoint');
endpt.x = temppt(1,1);
endpt.y = temppt(1,2);
width = abs(endpt.xstartpt.x)+0.00001;
height = abs(endpt.ystartpt.y)+0.00001;
hrect = findobj(h,'Tag','rangerect_fft');
set(hrect,'Position', [min(startpt.x, endpt.x), min(startpt.y, endpt.y), width, height]);
end
function BtnupFcn(h,evt)
set(h,'WindowButtonMotionFcn','');
set(h,'WindowButtonUpFcn','');
hrect = findobj(h,'Tag','rangerect_fft');
delete(hrect);
Filter_Analysis(startpt.x, endpt.x);
function Filter_Analysis(x1,x2)
lowfre = min(x1,x2);
highfre = max(x1,x2);
if
lowfre <= 0 || lowfre >= Fs/2
lowfre = 0.01;
end
if
highfre <= 0 || highfre >= Fs/2
highfre= Fs/2-0.01;
end
W1 = lowfre/(Fs/2);
W2 = highfre/(Fs/2);
filter_str = '(低通)';
switch filter_flag
case 1
[b,a] = butter(5,min(W1,W2)); % 低通，畫的框的左邊為截止頻率
filter_str = '(低通)';
case 2
[b,a] = butter(5,max(W1,W2),'high'); % 高通，畫的框的右邊為截止頻率
filter_str = '(高通)';
case 3
[b,a] = butter(5,[W1 W2]); % 帶通
filter_str = '(帶通)';
case 4
[b,a] = butter(5,[W1 W2],'stop'); % 帶阻
filter_str = '(帶阻)';
end
y = filter(b,a,yt);
figure
subplot(2,1,1)
plot(t,yt)
xlabel('時間t')
ylabel('信號y')
xlim([t(1) t(end)])
ylim([min(yt),max(yt)]);
title('原始信 號')
subplot(2,1,2)
plot(t,y)
xlabel('時間t')
ylabel('信號y')
xlim([t(1) t(end)])
ylim([min(y),max(y)]);
title(sprintf('濾波信號%s%.2fHz~%.2fHz',filter_str,lowfre,highfre))
end
end
end
end
subplot(1,3,[2,3]) % 頻率軸化為頻率
[X,Y] = meshgrid(t,5/(2*pi)./scale*Fs);
mesh(X,Y,abs(cw2))
view(0,90)
title('時頻圖')
xlabel('時間')
ylabel('頻率')
xlim([t(1) t(end)])
set(gca,'ylim',[0,max(max(Y))])
set(gca,'YScale','log')
set(gca,'YTick', [1:9,10:10:90,100:100:900,1000,2000])
function [Yf, f] = Spectrum_Calc(yt,Fs)
L = length(yt);
NFFT = 2^nextpow2(L);
Yf = fft(yt,NFFT)/L;
Yf = 2*abs(Yf(1:NFFT/2+1));
f = Fs/2*linspace(0,1,NFFT/2+1);
end
end
end
end
其中選框的代碼可以標準化，可以用來在用戶交互中用戶選擇范圍，代碼整理如下
set(gcf,'WindowButtonDownFcn',@BtndownFcn);
function BtndownFcn(h,evt)
temppt = get(gca,'CurrentPoint');
startpt.x = temppt(1,1);
startpt.y = temppt(1,2);
endpt = startpt;
height = 0.0001;
width = 0.0001;
rectangle('Position',[startpt.x, startpt.y, width, height],'Tag','rangerect'); set(h,'WindowButtonMotionFcn',@BtnmoveFcn);
set(h,'WindowButtonUpFcn',@BtnupFcn);
function BtnmoveFcn(h,evt)
temppt = get(gca,'CurrentPoint');
endpt.x = temppt(1,1);
endpt.y = temppt(1,2);
width = abs(endpt.x-startpt.x)+0.00001;
height = abs(endpt.y-startpt.y)+0.00001;
hrect = findobj('Tag','rangerect');
set(hrect,'Position',[min(startpt.x, endpt.x), min(startpt.y, endpt.y), width, height]);
end
function BtnupFcn(h,evt)
set(h,'WindowButtonMotionFcn','');
set(h,'WindowButtonUpFcn','');
hrect = findobj('Tag','rangerect');
delete(hrect);
% ProsessFcn(startpt,endpt);
% 選完框后對選擇的 范圍處理
end
end
