%function bandwidth
%對任意輸入信號求帶寬
% rect_a輸入信號
% dt_a輸入信號采樣周期
% threshold 帶寬門限 單位db
% ss_e輸出單邊功率譜密度信號
% f_high上限頻率
% f_low下限頻率
% bw_a輸入信號帶寬

function [ss_e,f_high,f_low,bw]=bandwidth(signal,dt,threshold)

%step1 evaluation of single sided esd

fs=1/dt;   %抽樣頻率
n=length(signal);  %抽樣點數，即fft長度
t=n*dt;          %時間窗口    
df=1/t;         %基本頻率
x=fft(signal);
x=x/n;        %matlab幅度譜轉為傅氏譜
ds_e=abs(x).^2/(df^2);     %雙邊功率譜
ss_e=2*ds_e(1:floor(n/2));  %單邊功率譜

%step2 evaluation of bandwidth

[epeak,index]=max(ss_e);
f_peak=index*df;
eth=epeak*10^(threshold/10);  %門限電平

%選擇上下限頻率
  imax=index;
  e0h=ss_e(index);
  while (e0h>eth)&(imax<=(n/2))
      imax=imax+1;
      e0h=ss_e(imax)
  end    %搜索上限頻率
  
  f_high=(imax-1)*df;
  
  imin=index;
  e0l=ss_e(index);
  while (e0l>eth)&(imin>1)&(index>1)
      imin=imin-1;
      e0l=ss_e(imin);
  end    %搜索下限頻率
  
  f_low=(min(index,imin)-1)*df;
  
  bw=f_high-f_low;   %信號帶寬
  
  fprintf('\nFrequency Bandwidth=%f  [Hz]\nHigh Frequency=%f [Hz]\nLow frequency=%f [Hz]',bw,f_high,f_low);
  
  %step3 graphical output
  
  figure(2)
    
   frequency=linspace(0,fs/2,length(ss_e));
   pf=plot(frequency,ss_e);
   set(pf,'linewidth',[2]);
   l1=line([f_high f_high],[min(ss_e),max(ss_e)]);
   set(l1,'color',[0 0 0],'linestyle',':')
   l1=line([f_low f_low],[min(ss_e) max(ss_e)]);
   set(l1,'color',[0 0 0],'linestyle',':')
   l1=line([f_low f_high],[eth eth]);
   set(l1,'linewidth',[2],'color','red','linestyle',':')
   axis([0.8*f_low 1.2*f_high -0.1*epeak 1.2*epeak]);
   ax=gca;
   set(ax,'fontsize',12);
   t=title('Frequency domain');
   set(t,'fontsize',14);
   x=xlabel('Frequency [Hz]');
   set(x,'fontsize',14)
   y=ylabel('Single_sided ESD  [V^2s/Hz]');
   set(y,'fontsize',14);