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
//
|
