MATLAB仿真程序
窗函數設計法：
Fs=100;
t=(1:100)/Fs;
s1=sin(2*pi*t*5);
s2=sin(2*pi*t*15);
s3=sin(2*pi*t*30);
s=s1+s2+s3;
figure(1);
subplot(2,2,1);
plot(t,s)%畫出信號的時域波形
xlabel('Time(seconds)')
ylabel('Time waveform')
title('原始信號的時域波形')%程序功能：畫出信號的頻譜圖。
S=fft(s,512);%對s進行快速傅立葉變換
w=(0:255)/256*(Fs/2);
subplot(2,2,2);
plot(w,abs(S(1:256)))%畫出信號的幅度圖
xlabel('Frequency(Hz)')
ylabel('幅度')
title('幅度譜')
axis([0 35 0 60]);
subplot(2,2,3);
plot(w,angle(S(1:256)))%畫出信號的相位圖
xlabel('Frequency(Hz)')
ylabel('相位')
title('相位譜')%程序功能：設計低通濾波器并畫出其頻譜圖：
%%%%%%%%%%%%%%低通
fb=10;
fc=13;%設置濾波器截止頻率
fs=100;
wb=2*pi*fb/fs;
ws=2*pi*fc/fs;
wc=0.5*(wb+ws);
tr_width=ws-wb;%過渡帶寬
%%%%%%%%%%%%%%高通
%fb=25;
%fc=22;%設置濾波器截止頻率
%fs=100;
%wb=2*pi*fb/fs;
%ws=2*pi*fc/fs;
%wc=0.5*(wb+ws);
%tr_width=wb-ws;%過渡帶寬
%%%%%%%%%%%%%%%%%%%帶通
%fc1=9;fb1=12;fb2=18;fc2=21;fs=100;
%wb1=2*pi*fb1/fs;ws1=2*pi*fc1/fs;wb2=2*pi*fb2/fs;
%ws2=2*pi*fc2/fs;wc1=0.5*(wb1+ws1);wc2=0.5*(wb2+ws2);
%tr_width=min((wb1-ws1),(ws2-wb2));
M=ceil(1.8*pi/tr_width);
n=[0:1:M-1];
hd=ideal_lp(wc,M);
%w_box=(boxcar(M))'
%w_box=(triang(M))'
%w_box=(hamming(M))'
w_box=(kaiser(M,2.5))'
h=hd.*w_box;
[H,w]=freqz(h,[1],1000,'whole');
H=(H(1:501))';
w=(w(1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(h,[1],w);
figure(2);
subplot(2,2,1);stem(n,hd);%理想脈沖響應
xlabel('n');
ylabel('hd(n)');
title('Ideal Impulse Response');
subplot(2,2,2);stem(n,w_box);%矩形窗
xlabel('n');
ylabel('w(n)');
title('Kaiser Window');
subplot(2,2,3);stem(n,h);%實際脈沖響應
xlabel('n');
ylabel('h(n)');
title('Actual Impulse Response');
subplot(2,2,4);plot(w*fs/(2*pi),db);%幅度響應（dB）
axis([0 40-50 0]);
xlabel('Frequency(Hz)');
ylabel('Decibels');
title('Magnitude Response in dB');
sf=filter(h,[1],s);%sf為濾濾波后的信號
figure(3);
subplot(2,2,1);
plot(t,sf)%畫出濾波后信號的時域波形
xlabel('Time(seconds)')
ylabel('Time waveform')
axis([0 1-1 1]);
title('濾波后信號的時域波形')
SF=fft(sf,512);%對sf進行快速傅里葉變換
w=(0:255)/256*(Fs/2);
subplot(2,2,2);
plot(w,abs(SF(1:256)))%畫出濾波后信號的幅度圖
xlabel('Frequency(Hz)')
ylabel('幅度譜')
title('濾波后信號的幅度譜');
subplot(2,2,3);
plot(w,angle(SF(1:256)))%畫出濾波后信號的相位圖
xlabel('Frequency(Hz)')
ylabel('相位譜')
title('濾波后信號的相位譜')%程序功能：對濾波前后信號進行比較
subplot(2,2,4);
plot(w,abs([S(1:256)'SF(1:256)']))%將濾波前后信號的幅度譜畫在一起
xlabel('Frequency(Hz)')
ylabel('Mag.of Fourier transform')
legend({'before','after'})%對兩個曲線進行區分命名
title('濾波前后信號對比')
頻率采樣設計法：
N=100;Fs=100; %數據總數和采樣頻率
fc1=12;
fc2=18;
n=[0:N-1];t=n/Fs;%時間序列
f1=5;f2=15;f3=30;
x=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);%輸入信號
%設計25階的低通濾波器，歸一化截止頻率為fc*2/Fs
b=fir1(25,0.2);
y1=fftfilt(b,x,256);%對數據進行濾波
n1=1:100;t1=t(n1);%選擇采樣點間隔
x1=x(n1);%與采樣點對應的輸入信號
subplot(2,2,1);plot(t1,x1);%繪制輸入信號
title('輸入信號');
S=fft(x,512);% 對x進行快速傅立葉變換
w=(0:255)/256*(Fs/2);
subplot(2,2,2);
plot(w,abs(S(1:256)))% 畫出信號的幅度圖
xlabel('Frequency (Hz)')
ylabel('幅度')
title('幅度譜')
n2=n1;t2=t(n2);%輸出信號，扣除了相位延遲N/2
y2=y1(n2);
subplot(2,2,3);plot(t2,y2);%繪制輸出信號
title('輸出信號');
xlabel('時間/s')
SF=fft(y1,512);% 對sf進行快速傅里葉變換
w=(0:255)/256*(Fs/2);
subplot(2,2,4);
plot(w,abs(SF(1:256)))% 畫出濾波后信號的幅度圖
xlabel('Frequency (Hz)')
ylabel('幅度譜')
title('濾波后信號的幅度譜');
figure(2);
[H,f]=freqz(b,1,512,Fs);
plot(f,20*log10(abs(H)));
DSP實現程序
FIR.c程序：
#include<math.h>
#define FIRNUMBER 25
#define PI 3.1415926
float InputWave();
float FIR();
FloatfHn[FIRNUMBER]={51455693,35039587,0,-43569197,-79886877,-91884544,-66619233,0,100778399,218060512,327780957,405671856,433829342,405671856,327780957,218060512,100778399,0,-66619233,-91884544,-79886877,-43569197,0,35039587,51455693};
float fXn[FIRNUMBER]={0.0};
float fInput,fOutput;
float fSignal1,fSignal2,fSignal3;
float fStepSignal1,fStepSignal2,fStepSignal3;
float f2PI;
int i;
float fIn[256],fOut[256];
int nIn,nOut;
main()
{
nIn=0;nOut=0;
f2PI=2*PI;
fSignal1=0.0;
fSignal2=0.0;
fStepSignal1=2*PI/20;
fStepSignal2=6*PI/20;
fStepSignal3=12*PI/20;
while(1)
{
fInput=InputWave();
fIn[nIn]=fInput;
nIn++;nIn%=256;
fOutput=FIR();
fOut[nOut]=-fOutput;
nOut++; /*break point*/
if(nOut>=256)
{
nOut=0;
}
}
}
float InputWave()
{
for(i=FIRNUMBER-1;i>0;i--)
fXn[i]=fXn[i-1];
fXn[0]=sin(fSignal1)+sin(fSignal2)+sin(fSignal3);
fSignal1+=fStepSignal1;
if(fSignal1>=f2PI) fSignal1-=f2PI;
fSignal2+=fStepSignal2;
if(fSignal2>=f2PI) fSignal2-=f2PI;
fSignal3+=fStepSignal3;
if(fSignal3>=f2PI) fSignal3-=f2PI;
return(fXn[0]);
}
float FIR()
{
float fSum;
fSum=0;
for(i=0;i<FIRNUMBER;i++)
{
fSum+=(fXn[i]*fHn[i]);
}
return(fSum);
}
Test.cmd程序：
-heap 100
-stack 400
-sysstack 400
-l rts55x.lib
MEMORY
{
PAGE 0:PROG(R):origin=0x80000,length=0x10000
PAGE 0:BOOT(R):origin=0x3FF000,length=0xFC0
PAGE 0:RESET(R):origin=0x3FFC0,length=0x2
//PAGE 0:VECTORS(R):origin=0x3FFC2,length=0x3E
PAGE 1:M0RAM(RW):origin=0x000000,length=0x400
PAGE 1:M1RAM(RW):origin=0x000400,length=0x400
PAGE 1:L0L1RAM(RW):origin=0x008000,length=0x2000
PAGE 1:H0RAM(RW):origin=0x3F8000,length=0x2000
}
SECTIONS
{
.reset:>RESET,PAGE=0
.pinit:>PROG,PAGE=0
.cinit:>PROG,PAGE=0
.text:>PROG,PAGE=0
.const:>L0L1RAM,PAGE=1
.bass:>L0L1RAM,PAGE=1
//
|
QQ好友和群
QQ空間
騰訊微博
騰訊朋友
收藏
淘帖
頂
踩
