clear all;                    	     				
clc;                              				
N=256;        %number of sample.                    			 	
fs=600;       %frequency of sample.                     			     
dt=1/fs;      %time interval of sample.                                  		
for k=1:N;    %gerenate integrated signal with f1,f2,f3.                
   f1=250;                         			 	
   f2=100;                         			 	
   f3=270;                       				
   y(k)=sin(2*pi*f1*k*dt)+sin(2*pi*f2*k*dt)+sin(2*pi*f3*k*dt);  
end

min(y);
max(y);
y=y*1024*8;
y=round(y);
n=1:256;
subplot(4,1,1)
plot(n,y,'b');
max(y);
A=y.';
%following code save data in vector into the file in the disc.
[filename,pathname] = uiputfile({'*.bin';'*.dat';'*.*'},'Save as'); 
fid = fopen([pathname,filename],'w'); 
if fid~=-1 %not equal to -1
  count=fprintf(fid,'%x\n',y+(y<0)*2^16);% write the data in vector 'y' to the file whose number is 'fid'. 
  if(fclose(fid)==-1) 
    error('Can not close the file!'); 
  end 
else 
  error('Can not open the file!'); 
end 

y
lp=200;       %frequency of cut-off.                            
wn1=2*lp/fs;  %歸一化的截止頻率。                            
[z1,p1,k1] = CHEBY1(3,0.5,wn1); %濾波器的零極點表示。          
[b1,a1] = CHEBY1(3,0.5,wn1);   %濾波器的傳輸函數表示。           
[H,f]=freqz(b1,a1,N/2,fs);     %濾波器的頻率響應
subplot(4,1,2);                %繪制濾波器的頻率響應。
plot(f,abs(H),'b'); 
grid;
%hold on;
%b1=[0.1598,0.4768,0.4768,0.1589];

%a1=[-0.1268,0.5239,-0.1257];
%[H1,f]=freqz(b2,a2,N/2,fs);
%plot(f,abs(H1),'r');

b1=b1/(8*1.0711)            %將參數按比列系數縮小。      		
a1=a1/(8*1.0711)          %為DSP程序做準備，MATLAB 程序不需要。           


yy1=filter(b1,a1,y);           %濾波: input signal--->y, output signal--->yy1;   	  

b1
a1

y=fft(y,N);                    %將原輸入信號做FFT變換。     	  
pyy=y.*conj(y);                %做原信號的256點的功率譜分析。"A.*A" 表示矩陣對應的元素相乘，而不矩陣相乘      	  
f=(0:(N/2-1));
for i=1:N/2-1;                 %還原為模擬頻率。    
   f(i)=f(i)*fs/N;
end
%figure(1);
subplot(4,1,3);
plot(f,pyy(1:N/2))             %繪制原輸入信號的功率譜圖。          
grid;


y=fft(yy1,N);                  %濾波后的輸出信號的FFT變換。    	 
pyy=y.*conj(y);                %計算輸出信號的功率譜密度。 
f=(0:(N/2-1));
for i=1:N/2-1;                 %還原為模擬頻率。    
   f(i)=f(i)*fs/N;
end
%figure(2);
subplot(4,1,4);
plot(f,pyy(1:N/2));           %繪制原輸出信號的功率譜圖。                     
grid;

%程序的整體思路：（1）：產生濾波的輸入信號。
%               （2）：根據截止頻率、波紋系數、通帶頻率設置濾波器。（[b1,a1] = CHEBY1(3,0.5,wn1);）   
%               （3）：對輸入信號濾波。（yy1=filter(b1,a1,y)）
%               （4）：對輸出信號做頻譜分析。（FFT）
%本實列要求使用MATLAB產生濾波前的3正弦合成信號，然后設計一個3階低通切比雪夫濾波器，
%濾除信號中的高頻信號，并使用頻譜分析方法，查看濾波效果。MATLAB仿真后，使用得到的濾波參數
%，進行DSP編程，在DSP中實現低通濾波，并使用CCS的頻譜分析功能，查看DSP的濾波效果。
%   由于本實列前兩個過程比較復雜，所以一般在MATLAB中驗證算法，并設計濾波的參數，再將此參數
%應用到DSP匯編語言程序中，從而簡化編程過程。然后將濾波后的信號通過DA輸出，使用示波器查看程序
%運行的正確性。也可以將濾波后的信號存儲在DSP的內部RAM中，使用CCS參看濾波后信號的時域和頻域圖。


%count=fprintf(fid,'%x\n',y+(y<0)*2^16);
%以上語句是把y向量中的數據以十六進制的格式寫到文件中去。
%對于y向量中的正數，就直接寫，如果是負數，對于32bit的要加上2^16后方能得到正確的結果。
%
%
%
%
%
%