clear all
Fs = 4000;            % Sampling frequency
Ts=1/Fs;    %sample period
Fi=100;    %message f or input signal
Wi=2*pi*Fi; %message angle velocity
n_fft=2^14;  %sample number 16384
t=(0:n_fft-1)*Ts; %sample total time
y=zeros(1,n_fft)
data=load(transpose('xxx.csv'));
y(1:n_fft)=data(2,1:n_fft); %data必須從非文字的那行(不一定是第2行)開始選取

subplot(2,1,1);     %在2X1中兩個圖繪製第一個圖

plot(t(1:100),y(1:100));

subplot(2,1,2);
%fft%
N=n_fft-1;
n=linspace(-N/2,N/2,n_fft); %-N/2~N/2分成n_fft個數
%Black man-Harris window function
w=0.35875+0.48829*cos(2*pi*n/(N-1))+0.14128*cos(4*pi*n/(N-1))+0.01168*cos(6*pi*n/(N-1));
y_w=y(1:n_fft).*w; %乘進去
Y_w=fft(y_w,n_fft);  %Discrete Fourier Transform

%plot spectrum%
Ys_w=2*abs(Y_w(1:n_fft/2));  % single side band
YM_w=max(Ys_w(1:n_fft/2));   % Ys_w message max magnitude
YL_w(1:n_fft/2)=20*log10(Ys_w(1:n_fft/2)/YM_w);  %magnitude of dB
f_w=Fs/2*linspace(0,1,n_fft/2); %f = Fs/2*(0~1分成n_fft/2個數)
%====第二張圖繪製====%
% subplot(2,1,2);  %在2X1中兩個圖繪製第二個圖
plot(f_w,YL_w); %f is the X axis and YL is the Y axis 
axis([0,Fs/2,-150,20]);     %determine the range of X and Y axis
xlabel('frequency(Hz)');   %the unit of X
ylabel('magnitude(dB)');   %the unit of Y
