;/*
=============================================================== =============*/
;/* Copyright (C 2004 YINXING TECHNOLOGY CO., LTD */
;/* All Rights Reserved. */
;/* ----------------------------------------------------------------------------*/
;/*
=============================================================== =============*/
;---------------------------------------------------------------
;
; USER Program Demo !
; It Start At 0x1800,and interrupt vector don't changed,
; At 200h !
;
; Please Don't modify PMST !
;
; Some data put in 080h DP=1
; SP may use system stack,so user needn't setup SP !
;
;------------------------------------------------------------------
.title "for test user program ... "
.mmregs
.global _c_int00
.ref fs_start,errno
op1valh .set 4140h ;floating point number 12.0
op1vall .set 0000h
op2valh .set 4140h ;floating point number 12.0
op2vall .set 0000h
initst0 .set 1800h ;set st0 initial number
initst1 .set 2900h ;set st1 initial number
.bss rlthm,1 ;result high mantissa,address is 80h .bss rltlm,1 ;result low mantissa,address is 81h
.bss rltsign,1 ;ressult sigh,address is 82h
.bss rltexp,1 ;result exponent,address is 83h
.bss op1hm,1 ;op1 high mantissa,address is 84h
.bss op1lm,1 ;op1 low mantissa,address is 85h
.bss op1se,1 ;op1 sigh and exp,address is 86h
.bss op2se,1 ;op2 sigh and exponent,address is 87h .bss op2hm,1 ;op2 high mantissa,address is 88h
.bss op2lm,1 ;op2 low mantissa,address is 89h
.bss op1_hsw,1 ;op1 packed high,address is 8ah .bss op1_lsw,1 ;op1 packed low,address is 8bh
.bss op2_hsw,1 ;op2 packed high,address is 8ch .bss op2_lsw,1 ;op2 packed low,address is 8dh
_c_int00:
stm #initst0,st0
stm #initst1,st1
rsbx C16
ld #op1lm,dp
ld #op1valh,a ;load float number
stl a,op1_hsw
ld #op1vall,a
stl a,op1_lsw
ld #op2valh,a ;load float number
stl a,op2_hsw
ld #op2vall,a
stl a,op2_lsw
nop
nop
nop ;1st breakpoint in CCS!
;----------- conversion of floating point format - unpack ------- dld op1_hsw,a ;load OP1 to acc a
sfta a,8
sfta a,-8
bc op1_zero,AEQ ;if op1 is 0,jump to special case
sth a,-7,op1se ;store sign and exponent to stack
stl a,op1lm ;store low mantissa
and #07Fh,16,a ;mask off sign and exp to get high mantissa add #080h,16,a ;add implied 1 to mantissa
sth a,op1hm ;store mantissa to stack
dld op2_hsw,a ;load OP2 to acc a
sfta a,8
sfta a,-8
bc op2_zero,AEQ ;if op2 is 0,jump to special case
sth a,-7,op2se ;store sign and exponent to stack
stl a,op2lm ;store low mantissa
and #07Fh,16,a ;mask off sign and exp to get high mantissa add #080h,16,a ;add implied 1 to mantissa
sth a,op2hm ;store mantissa to stack
nop
nop
nop ;2nd breakpoint in CCS !
;---------- judge the sign----------------
bitf op1se,#100h ;test the sign bit
bc testop2,NTC ;if is not negative jump to testop2
ld #0,a ;change the experssion to
dsub op1hm,a
dst a,op1hm ;store changed op1
testop2:
bitf op2se,#100h ;test the sign bit
bc compexp,NTC ;if is not negative jump to compexp
ld #0,a ;change the expression to
dsub op2hm,a
dst a,op2hm ;store changed op2
;--------- Exponent Comparison ------------
compexp:
ld op1se,a
and #00ffh,a ;mask off the sign bit
ld op2se,b
and #00ffh,b ;mask off the sign bit
sub a,b ;exp op2-exp op1 -> b
bc op1_gt_op2,BLT ;process op1 > op2
bc op2_gt_op1,BGT ;process op2 > op1
a_eq_b
dld op1hm,a
dadd op2hm,a ;add mantissa
bc res_zero,AEQ ;if result is zero process special case
ld op1se,b ;load exponent in preparation for normalizing normalize
sth a,rltsign ;Save signed mantissa on stack
abs a ;Create magnitude value of mantissa
sftl a,6 ;Pre–normalize adjustment of mantissa
exp a ;Get amount to adjust exp for normalization
norm a ;Normalize the result
st t,rltexp ;Store exp adjustment value
add #1,b ;Increment exp to account for implied carry
sub rltexp,b ;Adjust exponent to account for normalization normalized
stl b,rltexp ;Save result exponent on stack
bc underflow,BLEQ ;process underflow if occurs
sub #0ffh,b ;adjust to check for overflow
bc overflow,BGEQ ;process overflow if occurs
sftl a,-7 ;Shift right to place mantissa for splitting
stl a,rltlm ;Store low mantissa
and #07f00h,8,a ;Eliminate implied one
sth a,rlthm ;Save result mantissa on stack
;----------- Conversion of Floating Point Format- Pack --------- ld rltsign,9,a
and #100h,16,a ;Get the sign value
add rltexp,16,a ;Add the result exponent together
sftl a,7 ;shift the value to right place
dadd rlthm,a ;Add the result mantissa together
return_value
nop
nop
nop
nop ; 3th breakpoint in CCS !
b fs_start
op1_gt_op2
abs b ;if exp OP1 >= exp OP2 + 24 then return OP1
sub #24,b
bc return_op1,BGEQ
add #23,b ;restore exponent difference value
stl b,rltsign ;store exponent difference to be used as RPC
dld op2hm,a ;load OP2 mantissa
rpt rltsign ;normalize OP2 to match OP1
sfta a,-1
bd normalize ;delayed branch to normalize result
ld op1se,b ;load exponentvalue to prep for normalization dadd op1hm,a ;add OP1 to OP2
op2_gt_op1
sub #24,b ;if exp OP2 >= exp OP1 + 24 then return OP1 bc return_op2,BGEQ
add #23,b ;Restore exponent difference value
stl b,rltsign ;Store exponent difference to be used as RPC dld op1hm,a ;Load OP1 mantissa
rpt rltsign ;Normalize OP1 to match OP2
sfta a,-1
bd normalize ;Delayed branch to normalize result
ld op2se,b ;Load exponent value to prep for normalization dadd op2hm,a ;Add OP2 to OP1
op1_zero:
return_op2:
bd return_value
dld op2_hsw,a ;Put OP2 as result into A
op2_zero:
return_op1:
dld op1hm,a ;Load signed high mantissa of OP1
bc op1_pos,AGT ;If mantissa is negative .
neg a ;Negate it to make it a positive value
addm #100h,op1se ;Place the sign value back into op1_se op1_pos
sub #80h,16,a ;Eliminate implied one from mantissa
ld op1se,16,b ;Put OP1 back together in acc A as a result bd return_value
sftl b,7
add b,a
overflow
st #2,errno ;load error no
ld rltsign,16,a ;pack sign of result and #8000,16,a
or #0ffffh,a ;result low mantissa bd return_value
add #07f7fh,16,a ;result exponent underflow
st #1,errno ;load error no
b return_value
res_zero
bd return_value
sub a,a
nop
2.1 clc close all; n=0:15; p=8;q=2; x=exp(-(n-p.^2/q; figure(1; subplot(3,1,1; stem(n,x; title('exp(-(n-p^2/q,p=8,q=2'; xk1=fft(x,16; q=4; x=exp(-(n-p.^2/q; subplot(3,1,2; xk2=fft(x,16; stem(n,x; title('exp(-(n-p^2/q,p=8,q=4'; q=8; x=exp(-(n-p.^2/q;
xk3=fft(x,16; subplot(3,1,3; stem(n,x; title('exp(-(n-p^2/q,p=8,q=8';%时域特性figure(2; subplot(3,1,1; stem(n,abs(xk1; title('exp(-(n-p^2/q,p=8,q=2'; subplot(3,1,2; stem(n,abs(xk2; title('exp(-(n-p^2/q,p=8,q=4'; subplot(3,1,3; stem(n,abs(xk3; title('exp(-(n-p^2/q,p=8,q=8';%频域特性%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% p=8;q=8; figure(3; subplot(3,1,1; stem(n,x; title('exp(-(n-p^2/q,p=8,q=8';
xk1=fft(x,16; p=13; x=exp(-(n-p.^2/q; subplot(3,1,2; xk2=fft(x,16; stem(n,x; title('exp(-(n-p^2/q,p=13,q=8'; p=14; x=exp(-(n-p.^2/q; xk3=fft(x,16; subplot(3,1,3; stem(n,x; title('exp(-(n-p^2/q,p=14,q=8';%时域特性figure(4; subplot(3,1,1; stem(n,abs(xk1; title('exp(-(n-p^2/q,p=8,q=8'; subplot(3,1,2; stem(n,abs(xk2; title('exp(-(n-p^2/q,p=13,q=8'; subplot(3,1,3;
实验二 IIR数字滤波器设计及软件实现 1.实验目的 (1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法; (2)学会调用MATLAB信号处理工具箱中滤波器设计函数(或滤波器设计分析工具fdatool)设计各种IIR数字滤波器,学会根据滤波需求确定滤波器指标参数。 (3)掌握IIR数字滤波器的MATLAB实现方法。 (3)通过观察滤波器输入输出信号的时域波形及其频谱,建立数字滤波的概念。 2.实验原理 设计IIR数字滤波器一般采用间接法(脉冲响应不变法和双线性变换法),应用最广泛的是双线性变换法。基本设计过程是:①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标;②设计过渡模拟滤波器;③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。第六章介绍的滤波器设计函数butter、cheby1 、cheby2 和ellip可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。本实验要求读者调用如上函数直接设计IIR数字滤波器。 本实验的数字滤波器的MATLAB实现是指调用MATLAB信号处理工具箱函数filter对给定的输入信号x(n)进行滤波,得到滤波后的输出信号y(n)。 3. 实验内容及步骤 (1)调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st,该函数还会自动绘图显示st的时域波形和幅频特性曲线,如图1所示。由图可见,三路信号时域混叠无法在时域分离。但频域是分离的,所以可以通过滤波的方法在频域分离,这就是本实验的目的。 图1 三路调幅信号st的时域波形和幅频特性曲线 (2)要求将st中三路调幅信号分离,通过观察st的幅频特性曲线,分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器(低通滤波器、带通滤波器、高通滤波器)的通带截止频率和阻带截止频率。要求滤波器的通带最大衰减为0.1dB,阻带最小衰减为
3.6 常见的算法实现 在实际应用中虽然信号处理的方式多种多样,但其算法的基本要素却大多相同,在本节中介绍几种较为典型的算法实现,希望通过对这些例子(单精度,16bit )的分析,能够让大家熟悉DSP 编程中的一些技巧,在以后的工作中可以借鉴,达到举一反三的效果。 1. 函数的产生 在高级语言的编程中,如果要使用诸如正弦、余弦、对数等数学函数,都可以直接调用运行库中的函数来实现,而在DSP 编程中操作就不会这样简单了。虽然TI 公司提供的实时运行库中有一些数学函数,但它们所耗费的时间大多太长,而且对于大多数定点程序使用双精度浮点数的返回结果有点“大材小用”的感觉,因此需要编程人员根据自身的要求“定制”数学函数。实现数学函数的方法主要有查表法、迭代法和级数逼近法等,它们各有特点,适合于不同的应用。 查表法是最直接的一种方法,程序员可以根据运算的需要预先计算好所有可能出现的函数值,将这些结果编排成数据表,在使用时只需要根据输入查出表中对应的函数值即可。它的特点是速度快,但需要占用大量的存储空间,且灵活度低。当然,可以对上述查表法作些变通,仅仅将一些关键的函数值放置在表中,对任意一个输入,可根据和它最接近的数据采用插值方法来求得。这样占用的存储空间有所节约,但数值的准确度有所下降。 迭代法是一种非常有用的方法,在自适应信号处理中发挥着重要的作用。作为函数产生的一种方法,它利用了自变量取值临近的函数值之间存在的关系,如时间序列分析中的AR 、MA 、ARMA 等模型,刻画出了信号内部的特征。因为它只需要存储信号模型的参量和相关的状态变量,所以所占用的存储空间相对较少,运算时间也较短。但它存在一个致命的弱点,由于新的数值的产生利用了之前的函数值,所以它容易产生误差累积,适合精度要求不高的场合。 级数逼近法是用级数的方法在某一自变量取值范围内去逼近数学函数,而将自变量取值在此范围外的函数值利用一些数学关系,用该范围内的数值来表示。这种方法最大的优点是灵活度高,且不存在误差累积,数值精度由程序员完全控制。该方法的关键在于选择一个合适的自变量取值区间和寻找相应的系数。 下面通过正弦函数的实现,具体对上述三种方法作比较。 查表法较简单,只需要自制一张数据表,也可以利用C5400 DSP ROM 内的正弦函数表。 迭代法的关键是寻找函数值间的递推关系。假设函数采样时间间隔为T ,正弦函数的角频率为ω,那么可以如下推导: 令()()()T T ω?β?αω?-+=+sin sin sin 等式的左边展开为 T T side left ω?ω?sin cos cos sin _+= 等式的右边展开为 ()T T side right ω?βωα?sin cos cos sin _-+= 对比系数,可以得到1,cos 2-==βωαT 。令nT =?,便可以得到如下的递推式: [][][]21cos 2---=n s n s T n s ω
第一节 X=(-1)S×(1.M)×2E-127e=E-127 X=(-1)S×(1.M)×2E-1023 e=E-1023 我承认以前对这俩公式避之不及不予深究努力自己说服自己而未能得逞,部分原因是跟“移码与真值的关系”扯上关系,这“移码与真值的关系”想搞清先得把引入移码的充分理由给我个说法,不幸玩过头正事误了。上回说了“补码省心移码悦目”能算是今时不同往日了吧,现在轮到对IEEE754浮点数规格化表示法杀无赦去死吧。 首先,“IEEE规格化形式”是对“传统规格化形式”进一步严格要求来的。 IEEE规格化形式唯一,而浮点数记法多种多样。 (1.75)10=1.11×20 (IEEE规格化表示)=0.111×21 (传统规格化表示) =0.0111×22=0.00111×23 其次,既然IEEE想到对“传统规格化形式”进一步修订当然有目的,你以为作无用功呐,关键目的是什么? 规格化的目的同理。修改阶码同时左右移小数点使尾数域最高有效位固定为1,尾数就以ta所可能变化成的最大形式出现,即使遭遇类似截断的操作仍可保持尽可能高的精度。 有类错误我这种大秀逗极善于犯!就是不理会左右关系不经过大脑直接作问题少女状问很白的问题:“‘移码和真值的关系’是E=27(或210)+X,那X=E-27(或210),在怎么着里面数该是128(或1024),咋是127(或1023)?” 当E=M=全0 E(移码)=全0,对应真值-128 M(补码)=全0,对应真值0 E=M=全0,真值X=0-128=0 结合符号位S 为0或1分正零和负零 当E=全1,M=全0 E(移码)=全1,对应真值+127 M(补码)=全0,对应真值0 E=全1,M=全0,真值X=0127=∞ 结合符号位S 为0或1分+∞和-∞ 要除去表示零和无穷大这2种特殊情况 指数偏移值不选128(10000000),而选127(01111111) 对IEEE32位规格化浮点数 8位移码(隐含1位符号位)原本表示范围是-128 →+127 (除去全1(+127)全0(-128)剩下-127 →+126 ???) 实际可用指数值(即阶码真值)e范围是-126→+127 加上偏移值后,阶码E的范围变为1→254 以10的幂表示,绝对值的范围是10-38→1038 假设由S,E,M三个域组成的一个32位二进制字所表示的非零规格化浮点数x,真值表示为:x=(-1)s×(1.M)×2E-128 它所表示的规格化的最大正数、最小正数、最大负数、最小负数是多少? 第二节 1、什么是IEEE754标准 用来规范化浮点数,其格式是
数字信号处理的步骤与注意事项,并编写1024个采样点的FFT C语言程序1. 数字信号处理 1.1 数字信号处理概述 数字信号处理是研究如何用数字或符号序列来表示信号以及如何对这些序列进行处理的一门学科。随着计算机技术的发展,数字信号处理技术得到了越来越广泛的应用,它已成为现代科学技术必不可少的工具。数字信号是数据序列,其处理实际上就是进行各种数学运算,如加、减、乘以及各种逻辑运算等等。因此,数字信号处理既可以是软件处理也可以是硬件处理。所谓软件处理,就是对所需要的运算编制程序,然后在计算机上实现,其处理灵活、方便。所谓硬件处理,就是用加法器、乘法器、延时器以及它们的各种组合来构成数字电路,以实现所需要的运算。硬件处理显然不如软件处理灵活方便,但能对数字信号进行实时处理。近年来日益广泛采用的各种数字信号处理器(如TI TMS320系列、Philps Trimedia系列等)可以认为是软硬件处理方式的结合,这种处理时用数字信号处理芯片以及存储器来组成硬件电路,所需要的运算靠特定的汇编语言编程来实现。因此,采用数字信号处理器既方便灵活,又能做到实时处理,所以数字信号处理器(DSP)已经越来越广泛地应用于包括通信在内的各个领域之中。 1.2 数字信号处理的优点 (1)精度高 数字系统的特性不因环境的变化而变化,计算精度是模拟系统所无法相比的,运算位数由8位提高到16位、32位、64位。 (2)可靠性高 模拟系统中各种参数受温度、环境影响较大,因而易出现感应、杂散效应,甚至会出现震荡等等;而数字系统受温度、环境影响较小。模拟信号受到干扰即产生失真,而数字信号由于只有两种状态,因此,所受的干扰只要在一定范围以内,就不会产生影响,这就是说,数字信号抗干扰能力强。另外,如果用数字信号进行传输,在中继站还可以再生。总的说来,信号的数字处理可靠性高。(3)灵活性强 可以通过改变数字信号系统的参数来改变系统的性能。数字信号的灵活性还表现在可以利用一套计算设备同时处理多路相互独立的信号,即所谓的“时分复用”,这在数字电话系统中是非常有用的技术。 (4)便于大规模集成化 数字部件具有高度的规范性,易于实现大规模集成化和大规模生产,数字系统体积小、重量轻。 (5)数字信号便于加密处理 由于数字信号实际上为数据序列,因此便于加密运算处理。
基于DSP汇编语言的高精度三角函数仿真及实现 【摘要】在现今DSP软件工程设计中,广泛采用高级语言(如C语言)直接调用三角函数进行计算。然而,汇编语言的稳定性、可读性和执行效率均优于高级语言,因缺少可供直接调用的三角函数库,其开发周期长,且计算精度难以保证,在工程中推广较少。本文就基于DSP汇编语言,提出一种高精度三角函数算法。 【关键词】三角函数传统查表法线性查表法MATLAB 仿真 一、引言 随着现代科技的高速发展,带三角函数的大型DSP复杂算法被广泛应用于各个领域,高级语言(如C语言)编程实现时,其稳定性、可读性、执行效率均不及汇编语言。然而,汇编算法没有标准的三角函数库,需软件开发人员自己编写,这就延长了软件的开发周期,同时,其计算精度难以得到保证。本文就基于DSP汇编语言,提出一种占用资源低且精度较高的线性查表法,实现三角函数,并通过MATLAB 仿真与传统查表法进行对比分析。 二、查找表建模 假设以等步长step定义了点数为N的三角函数查找表,
表中地址n处的值为(x0,y0),若要计算某x处的三角函数值y,用传统查找表实现,即近似于地址n处的查表值y0。其查找生成方式为:TAB_X(n)=Fx(n*step),n=0,1,2……N-1。其中TAB_X为三角函数表,Fx代表相应的标准三角函数。 而用线性查找表实现,即近似于(x0,y0)点切线上x 处所对应的值y’。设切线的斜率为C,则用公式可表示为:y≈y’=y0+C×=y0+C*(x-x0) 由公式可知,线性查找表需要包含x0处的三角函数值y0和切线斜率C,因此,采取y0值与C值交替存放方式,即为:TAB_X(n)=Fx(n*step),n=0,2,4……N-2;TAB_X (n)=Fc(n*step),n=1,3,5……N-1。其中Fc代表斜率函数,其对应关系见表1。 三、工程实现方法 对正弦和余弦函数,其函数本身与切线斜率正好互补,因此可共用一张查找表。此外,因正弦和余弦函数互余,可通过或将其定义域映射到0~/4。对反三角函数,其计算的值域范围为0~/2,随着x0值的增加,其斜率值单调递增,这将会导致图1中的y’值逐渐偏离y值,增大误差。因此,采取将反三角函数/4~/2的值域映射到/4~0的值域中进行 计算。其中,反正弦、反余弦、反正切函数分别通过、、进行映射,这样,反正弦和反余弦函数的定义域范围为0~
数字信号处理实验实验一离散时间信号与系统及MA TLAB实现 1.单位冲激信号: n = -5:5; x = (n==0); subplot(122); stem(n, x); 2.单位阶跃信号: x=zeros(1,11); n0=0; n1=-5; n2=5; n = n1:n2; x(:,n+6) = ((n-n0)>=0); stem(n,x); 3.正弦序列: n = 0:1/3200:1/100; x=3*sin(200*pi*n+1.2); stem(n,x); 4.指数序列 n = 0:1/2:10; x1= 3*(0.7.^n); x2=3*exp((0.7+j*314)*n); subplot(221); stem(n,x1); subplot(222); stem(n,x2); 5.信号延迟 n=0:20; Y1=sin(100*n); Y2=sin(100*(n-3)); subplot(221); stem(n,Y1); subplot(222); stem(n,Y2);
6.信号相加 X1=[2 0.5 0.9 1 0 0 0 0]; X2=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7]; X=X1+X2; stem(X); 7.信号翻转 X1=[2 0.5 0.9 1]; n=1:4; X2=X1(5-n); subplot(221); stem(n,X1); subplot(222); stem(n,X2); 8.用MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。a=[-2 0 1 -1 3]; b=[1 2 0 -1]; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel('n'); ylabel('幅度'); 9.用MA TLAB计算差分方程 当输入序列为时的输出结果。 N=41; a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)]; k=0:1:N-1; y=filter(a,b,x); stem(k,y) xlabel('n'); ylabel('幅度') 10.冲激响应impz N=64; a=[0.8 -0.44 0.36 0.22];
一、实验目的 1. 通过本次实验回忆并熟悉MATLAB这个软件。 2. 通过本次实验学会如何利用MATLAB进行序列的简单运算。 3. 通过本次实验深刻理解理论课上的数字信号处理的一个常见方法——对时刻n的样本附近的一些样本求平均,产生所需的输出信号。 3. 通过振幅调制信号的产生来理解载波信号与调制信号之间的关系。 二、实验内容 1. 编写程序在MATLAB中实现从被加性噪声污染的信号中移除噪声的算法,本次试验采用三点滑动平均算法,可直接输入程序P1.5。 2. 通过运行程序得出的结果回答习题Q1.31-Q1.33的问题,加深对算法思想的理解。 3. 编写程序在MATLAB中实现振幅调制信号产生的算法,可直接输入程序P1.6。 4. 通过运行程序得出的结果回答习题Q1.34-Q1.35的问题,加深对算法思想的理解。 三、主要算法与程序 1. 三点滑动平均算法的核心程序: %程序P1.5 %通过平均的信号平滑 clf; R=51; d=0.8*(rand(R,1)-0.5);%产生随噪声 m=0:R-1; s=2*m.*(0.9.^m);%产生为污染的信号 x=s+d';%产生被噪音污染的信号 subplot(2,1,1); plot(m,d','r-',m,s,'g--',m,x,'b-.');
xlabel('时间序号n');ylabel('振幅'); legend('d[n]','s[n]','x[n]'); x1=[0 0 x];x2=[0 x 0];x3=[x 0 0]; y=(x1+x2+x3)/3; subplot(2,1,2); plot(m,y(2:R+1),'r-',m,s,'g--'); legend('y[n]','s[n]'); xlabel('时间序号n');ylabel('振幅'); 2. 振幅调制信号的产生核心程序:(由于要几个结果,因此利用subplot函数画图) %程序P1.6 %振幅调制信号的产生 n=0:100; m=0.1;fH=0.1;fL=0.01; m1=0.3;fH1=0.3;fL1=0.03; xH=sin(2*pi*fH*n); xL=sin(2*pi*fL*n); y=(1+m*xL).*xH; xH1=sin(2*pi*fH1*n); xL1=sin(2*pi*fL1*n); y1=(1+m1*xL).*xH; y2=(1+m*xL).*xH1; y3=(1+m*xL1).*xH; subplot(2,2,1); stem(n,y); grid; xlabel('时间序号n');ylabel('振幅');title('m=0.1;fH=0.1;fL=0.01;'); subplot(2,2,2); stem(n,y1); grid; xlabel('时间序号n');ylabel('振幅');title('m=0.3;fH=0.1;fL=0.01;'); subplot(2,2,3); stem(n,y2); grid; xlabel('时间序号n');ylabel('振幅');title('m=0.3;fH=0.3;fL=0.01;'); subplot(2,2,4); stem(n,y3); grid;
摘要 泰勒级数展开法作为一种数学方法,在科研和平时的数据处理 方面应用的很广泛。尤其是在通信、仪器仪表和工业控制等领域应 用更为广泛。在科技高速发展的今天,对函数的计算不仅要求有很 高的精度,还对计算的时间又很高的要求,必须在很短的时间内完 成数据的处理,否则根本不能完成大批量数据的实时性计算和处理。介于DSP芯片运算速度快的特点,用DSP芯片完成这些算法已越来 越受到重视。 产生正弦波,分别是查表法和泰勒级数展开法。查表法应用于 精度要求不很高的场合,而泰勒级数展开法是一种比查表法更为有 效的方法。它能精确的计算出一个角度的正弦和余弦值,且占用的 储存空间较小,体现了它的优越性。 关键词: DSP 泰勒级数正弦波
目录 摘要 .....................................................................................................I 1、正余弦信号发生器的实现原理 .. (1) 1.1、正弦波信号发生器 (1) 2、正弦波的实现 (2) 2.1、计算一个角度的正弦值 (2) 2.2、计算一个角度的余弦值 (4) 3、正弦波的实现 (8) 4、链接文件 (10) 5、调试结果 (12) 总结 (13) 参考文献 (14)
1、正余弦信号发生器的实现原理 1.1、正弦波信号发生器 泰勒级数展开法是根据泰勒展开式进行计算来实现正弦信号,它能精确地计算出一个角度的正弦和余弦值,且只需要较小的存储空间。 本次主要用泰勒级数展开法来实现正弦波信号。 正弦函数和余弦函数可以展开成泰勒级数,其表达式: 3579 sin()3!5!7!9!x x x x x x =-+-+- 2468 cos()12!4!6!8! x x x x x =-+-+- 取泰勒级数的前5项,得近似计算式: 3579 2222 sin()3!5!7!9! 111123456789(((()))) x x x x x x x x x x x =-+-+ =----???? 2468 2222 cos()12!4!6!8! 11112345678 ((())) x x x x x x x x x =-+-+ =----??? 递推公式: sin(nx ) = 2cos(x )sin[(n -1)x ]-sin[(n -2)x ] cos(nx ) = 2cos(x )sin[(n -1)x ]-cos[(n -2)x ] 由递推公式可以看出,在计算正弦和余弦值时,需要已知cos(x )、sin(n -1)x 、sin(n -2)x 和cos(n -2)x 。
数字信号处理课程设计 设计题目: 姓名: 学号: 院系班级: 组次: 指导教师: 时间:2015年11月21日——2015年12月6日
摘要 基于 MATLAB 的图像边缘检测算法的研究和实现 图像边缘是图像的最基本的特征。所谓边缘,就是指图像局部强度变化最明显的部分,存在于区域与区域、目标与目标、目标与背景、基元与基元之间,包含有图像处理中用于识别的关键信息。边缘检测是数字图像处理中,最基础也是最重要的环节之一。本文介绍了六种经典的边缘检测算子,包括 Roberts 算子,Sobel 算子,Canny算子,Prewitt 算子,LOG 算法。并且利用 MATLAB 系统所提供的相关函数等,对同一副图像结合用这些不同的算子分别进行处理,分析并得到他们处理图像的特点。比较传统的边缘检测算子,因为是基于图像函数的一阶导数进行考察的,因而它们具有共同的特点是计算简单、速度较快,但是对噪声都比较敏感。LOG 算法和 Canny算法,都是先对图像进行平滑去噪,抗噪性能较好,但是会损失一些边缘信息,其中 LOG算法比较适合处理渐变灰度图像,而 Canny 算子更适合处理阶跃型边缘图像。小波变换边缘检测法,则能够很好的保留图像的边缘信息,更适合处理小阵列图像。 关键词: MATLAB;图像处理;边缘检测;微分算子
目录 第一章绪论 (4) 1.1设计目的与要求 (4) 1.2叙述国内外研究动态 (5) 第二章软件设计- 基于MatLab的边缘检测算法 (6) 2.1 MatLab简介 (6) 2.2边缘检测算法原理 (7) 2.2.1 Roberts 边缘算子 (7) 2.2.2 Sobel 边缘算子 (8) 2.2.3 Prewitt 边缘算子 (8) 2.2.4 Log 边缘算子 (8) 2.2.5 Canny 边缘算子 (8) 2.3边缘检测算法--测试程序 (9) 第三章实验结果及分析 (13) 3.1 Roberts算子检测图像边缘的实现 (13) 3.2 Sobel算子检测图像边缘的实现 (14) 3.3 Prewitt算子检测图像边缘的实现 (15) 3.4高斯一拉普拉斯LOG算子检测图像边缘的实现 (16) 3.5 Canny算子检测图像边缘的实现 (17) 第四章总结与心得体会 (20) 参考文献 (21) 致谢 (22)
实验一 MATLAB 仿真软件的基本操作命令和使用方法 实验容 1、帮助命令 使用 help 命令,查找 sqrt (开方)函数的使用方法; 2、MATLAB 命令窗口 (1)在MATLAB 命令窗口直接输入命令行计算3 1)5.0sin(21+=πy 的值; (2)求多项式 p(x) = x3 + 2x+ 4的根; 3、矩阵运算 (1)矩阵的乘法 已知 A=[1 2;3 4], B=[5 5;7 8],求 A^2*B
(2)矩阵的行列式 已知A=[1 2 3;4 5 6;7 8 9],求A (3)矩阵的转置及共轭转置 已知A=[1 2 3;4 5 6;7 8 9],求A' 已知B=[5+i,2-i,1;6*i,4,9-i], 求B.' , B' (4)特征值、特征向量、特征多项式 已知A=[1.2 3 5 0.9;5 1.7 5 6;3 9 0 1;1 2 3 4] ,求矩阵A的特征值、特征向量、特征多项式;
(5)使用冒号选出指定元素 已知:A=[1 2 3;4 5 6;7 8 9];求A 中第3 列前2 个元素;A 中所有列第2,3 行的元素; 4、Matlab 基本编程方法 (1)编写命令文件:计算1+2+…+n<2000 时的最大n 值;
(2)编写函数文件:分别用for 和while 循环结构编写程序,求 2 的0 到15 次幂的和。
5、MATLAB基本绘图命令 (1)绘制余弦曲线 y=cos(t),t∈[0,2π]
(2)在同一坐标系中绘制余弦曲线 y=cos(t-0.25)和正弦曲线 y=sin(t-0.5), t∈[0,2π] (3)绘制[0,4π]区间上的 x1=10sint 曲线,并要求: (a)线形为点划线、颜色为红色、数据点标记为加号; (b)坐标轴控制:显示围、刻度线、比例、网络线 (c)标注控制:坐标轴名称、标题、相应文本; >> clear;
太原理工大学DSP原理及应用 课程设计报告 专业班级:通信0802 姓名:邢剑卿 学号:2008001316
2010年12月30日 一、设计目的 学会用ccs 集成开发软件,在开发环境下完成工程项目创建,程序编辑,编译,链接,调试和数据分析。 二、设计内容 编写程序,利用ccs 软件产生正弦波 三、设计原理 正弦波信号发生器已被广泛地应用于通信、仪器仪表和工业控制等领域的信号处理系统中。 通常有两种方法可以产生正弦波,分别为查表法和泰勒级数展开法。 查表法是通过查表的方式来实现正弦波,主要用于对精度要求不很高的场合。 泰勒级数展开法是根据泰勒展开式进行计算来实现正弦信号,它能精确地计算出一个角度的正弦和余弦值,且只需要较小的存储空间。 本次主要用泰勒级数展开法来实现正弦波信号。 产生正弦波的算法 正弦函数和余弦函数可以展开成泰勒级数,其表达式: 取泰勒级数的前5项,得近似计算式: 递推公式: sin(nx ) = 2cos(x )sin[(n -1)x ]-sin[(n -2)x ] cos(nx ) = 2cos(x )sin[(n -1)x ]-cos[(n -2)x ] 由递推公式可以看出,在计算正弦和余弦值时,需要已知cos(x )、sin(n -1)x 、sin(n -2)x 和cos(n -2)x 。 -+-+-=!9!7!5!3)sin(9753x x x x x x -+-+-=!8!6!4!21)cos(8642x x x x x ))))((((981761541321 ! 9!7!5!3)sin(2 2229753?-?-?-?-=+-+-=x x x x x x x x x x x ))) (((87165143121 ! 8!6!4!21)cos(2 2228642?-?-?--=+-+-=x x x x x x x x x
浮点数的表示和基本运算 1 浮点数的表示 通常,我们可以用下面的格式来表示浮点数 S P M 其中S是符号位,P是阶码,M是尾数 对于IBM-PC而言,单精度浮点数是32位(即4字节)的,双精度浮点数是64位(即8字节)的。两者的S,P,M所占的位数以及表示方法由下表可知 S P M表示公式偏移量 1823(-1)S*2(P-127)*1.M127 11152(-1)S*2(P-1023)*1.M1023 以单精度浮点数为例,可以得到其二进制的表示格式如下 S(第31位)P(30位到 23位) M(22位到 0位) 其中S是符号位,只有0和1,分别表示正负;P是阶码,通常使用移码表示(移码和补码只有符号位相反,其余都一样。对于正数而言,原码,反码和补码都一样;对于负数而言,补码就是其绝对值的原码全部取反,然后加1.) 为了简单起见,本文都只讨论单精度浮点数,双精度浮点数也是用一样的方式存储和表示的。 2 浮点数的表示约定 单精度浮点数和双精度浮点数都是用IEEE754标准定义的,其中有一些特殊约定。 (1) 当P = 0, M = 0时,表示0。 (2) 当P = 255, M = 0时,表示无穷大,用符号位来确定是正无穷大还是负无穷大。
(3) 当P = 255, M != 0时,表示NaN(Not a Number,不是一个数)。 当我们使用.Net Framework的时候,我们通常会用到下面三个常量 Console.WriteLine(float.MaxValue); // 3.402823E+38 Console.WriteLine(float.MinValue); //-3.402823E+38 Console.WriteLine(float.Epsilon); // 1.401298E-45 //如果我们把它们转换成双精度类型,它们的值如下 Console.WriteLine(Convert.ToDouble(float.MaxValue)); // 3.40282346638529E+38 Console.WriteLine(Convert.ToDouble(float.MinValue)); //-3.40282346638529E+38 Console.WriteLine(Convert.ToDouble(float.Epsilon)); // 1.40129846432482E-45 那么这些值是如何求出来的呢? 根据上面的约定,我们可以知道阶码P的最大值是11111110(这个值是254,因为255用于特殊的约定,那么对于可以精确表示的数来说,254就是最大的阶码了)。尾数的最大值是11111111111111111111111。 那么这个最大值就是:0 11111110 11111111111111111111111。 也就是 2(254-127) * (1.11111111111111111111111)2 = 2127 * (1+1-2-23) = 3.40282346638529E+38 从上面的双精度表示可以看出,两者是一致的。最小的数自然就是- 3.40282346638529E+38。 对于最接近于0的数,根据IEEE754的约定,为了扩大对0值附近数据的表示能力,取阶码P = -126,尾数 M = (0.00000000000000000000001)2 。此时该数的二进制表示为:0 00000000 00000000000000000000001 也就是2-126 * 2-23 = 2-149 = 1.40129846432482E-45。这个数字和上面的Epsilon 是一致的。 如果我们要精确表示最接近于0的数字,它应该是 0 00000001 00000000000000000000000 也就是:2-126 * (1+0) = 1.17549435082229E-38。 3 浮点数的精度问题 浮点数以有限的32bit长度来反映无限的实数集合,因此大多数情况下都是一个近似值。同时,对于浮点数的运算还同时伴有误差扩散现象。特定精度下看似
1. 矩形窗: 程序代码: wp=0.2*pi; wst=0.3*pi; tr_width=wst-wp; N(1)=ceil(1.8*pi/tr_width)+1; w_boxcar=boxcar(N(1))'; N(2)=ceil(6.2*pi/tr_width)+1; w_hanning=hanning(N(2))'; N(3)=ceil(6.6*pi/tr_width)+1; w_hamming=hamming(N(3))'; N(4)=ceil(11*pi/tr_width)+1; w_blackman=blackman(N(4))'; N(5)=ceil((50-7.95)/(2.285*tr_width)+1); w_kaiser=kaiser(N(5),0.1102*(50-8.7))'; n=0:(N(1)-1); wc=(wp+wst)/2; alpha=(N(1)-1)/2; hd=(wc/pi)*sinc(wc/pi*(n-alpha)); h=hd.*w_boxcar; figure(1); subplot(221);stem(n,hd,'filled'); axis tight;xlabel('n');ylabel('hd(n)'); [Hr,w1]=zerophase(h); subplot(222);plot(w1/pi,Hr); axis tight;xlabel('\omega/\pi');ylabel('H(\omega)'); subplot(223);stem(n,h,'filled'); axis tight;xlabel('n');ylabel('h(n)'); [H,w]=freqz(h,1); subplot(224);plot(w/pi,20*log10(abs(H)/max(abs(H)))); xlabel('\omega/\pi');ylabel('db'); grid on 程序结果:
单双精度浮点数的IEEE标准格式 目前大多数高级语言(包括C)都按照IEEE-754标准来规定浮点数的存储格式,IEEE754规定,单精度浮点数用4字节存储,双精度浮点数用 8字节存储,分为三个部分:符号位、阶和尾数。阶即指数,尾数即有效小数位数。单精度格式阶占8位,尾数占24位,符号位1位,双精度则为11为阶,53 位尾数和1位符号位,如下图所示: 31 30 23 22 0 63 62 52 51 0 细心的人会发现,单双精度各部分所占字节数量比实际存储格式都了一位,的确是这样,事实是,尾数部分包括了一位隐藏位,允许只存储23位就可以表示24位尾数,默认的1位是规格化浮点数的第一位,当规格化一个浮点数时,总是调整它使其值大于等于1而小于2,亦即个位总是为1。例如1100B,对其规格化的结果为1.1乘以2的三次方,但个位1并不存储在23位尾数部分内,这个1是默认位。 阶以移码的形式存储。对于单精度浮点数,偏移量为127(7FH),而双精度的偏移量为1023(3FFH)。存储浮点数的阶码之前,偏移量要先加到阶码上。前面例子中,阶为2的三次方,在单精度浮点数中,移码后的结果为127+3即130(82H),双精度为1026(402H)。 浮点数有两个例外。数0.0存储为全零。无限大数的阶码存储为全1,尾数部分全零。符号位指示正无穷或者负无穷。 下面举几个例子:
所有字节在内存中的排列顺序,intel的cpu按little endian顺序,motorola 的cpu按big endian顺序排列。
IEEE754标准的一个规格化 32位浮点数x的真值可表示为 x=(-1)^S*(1.M)*2^(E-127)e=E-127 31 30 23 0 |S | E |M | [例1]若浮点数x的754标准存储格式为(41360000)16,求其浮点数的十进制数值。 解:将16进制展开后,可得二进制数格式为 0 100,0001,0 011,0110,0000,0000,0000,0000 S E M 指数e=100,0001,0-01111111=00000011=(3)10 包含隐藏位1的尾数1.M=1.011,0110,0000,0000,0000,0000 于是有x=(-1)^0*(1.M)*2^(E-127) =+(1.011011)2*2^3 =(11.375)10 [例2]将数(20.59375)10转化为754标准的32位浮点数的二进制存储格式。解:首先分别将整数部分和小数部分转换成二进制 (20.59375)10=+(10100.10011)2 然后移动小数点使其在1,2位之间 10100.10011=1.010010011*2^4 e=4 于是得到:S=0,E=e+127=131,M=010010011 最后得到32位浮点数的二进制存储格式为 0 100,0001,1 010,0100,1100,0000,0000,0000 =(41A4C000)16 从存储结构和算法上来讲,double和float是一样的,不一样的地方仅仅是float是32位的,double是64位的,所以double能存储更高的精度。 任何数据在内存中都是以二进制(0或1)顺序存储的,每一个1或0被称为1位,而在 x86CPU上一个字节是8位。比如一个16位(2字节)的 short int型变量的值是1000,那么它的二进制表达就是:00000011 11101000。由于Intel CPU的架构原因,它是按字节倒序存储的,那么就因该是这样:11101000 00000011,这就是定点数1000在内存中的结构。 目前C/C++编译器标准都遵照IEEE制定的浮点数表示法来进行float,double运算。这种结构是一种科学计数法,用符号、指数和尾数来表示,底数定为2——即把一个浮点数表示为尾数乘以2的指数次方再添上符号。下面是具体的规格: ````````符号位阶码尾数长度 float 1 8 23 32 double 1 11 52 64
实验一熟悉MATLAB环境 一、实验目的 (1)熟悉MATLAB的主要操作命令。 (2)学会简单的矩阵输入和数据读写。 (3)掌握简单的绘图命令。 (4)用MATLAB编程并学会创建函数。 (5)观察离散系统的频率响应。 二、实验内容 认真阅读本章附录,在MATLAB环境下重新做一遍附录中的例子,体会各条命令的含义。在熟悉了MATLAB基本命令的基础上,完成以下实验。 上机实验内容: (1)数组的加、减、乘、除和乘方运算。输入A=[1 2 3 4],B=[3 4 5 6],求C=A+B,D=A-B,E=A.*B,F=A./B,G=A.^B并用stem语句画出A、B、C、D、E、F、G。 实验程序: A=[1 2 3 4]; B=[3 4 5 6]; n=1:4; C=A+B;D=A-B;E=A.*B;F=A./B;G=A.^B; subplot(4,2,1);stem(n,A,'fill');xlabel ('时间序列n');ylabel('A'); subplot(4,2,2);stem(n,B,'fill');xlabel ('时间序列n ');ylabel('B'); subplot(4,2,3);stem(n,C,'fill');xlabel ('时间序列n ');ylabel('A+B'); subplot(4,2,4);stem(n,D,'fill');xlabel ('时间序列n ');ylabel('A-B'); subplot(4,2,5);stem(n,E,'fill');xlabel ('时间序列n ');ylabel('A.*B'); subplot(4,2,6);stem(n,F,'fill');xlabel ('时间序列n ');ylabel('A./B'); subplot(4,2,7);stem(n,G,'fill');xlabel ('时间序列n ');ylabel('A.^B'); 运行结果:
目录一.引言 二.实验原理 三.FFT基本结构 (1)信号流图 (2)软件程序流图 四.实验程序 五.调试过程与步骤 六.实验结果 七.结果分析 八.遇到的问题及解决办法 九.实验体会
实验题目:快速傅立叶变换(FFT)算法一.引言 众所周知,FFT 是离散傅立叶变换(DFT)的一种快速算法。由于计算DFT 时一次复数乘法需用四次实数乘法和二次实数加法;一次复数加法则需二次实数加法。每运算一个X(k)需要4N 次复数乘法及2N+2(N-1)=2(2N-1)次实数加法。所以整个DFT运算总共需要4N^2 次实数乘法和N*2(2N-1)=2N(2N-1)次实数加法。如此一来,计算时乘法次数和加法次数都是和N^2 成正比的,当N 很大时,运算量是可观的,因而需要改进对DFT 的算法提高运算速度。 我们观察DFT的系数特性,可发现其对称性和周期性,我们可以将DFT 运算中有些项合并。因此,FFT就孕育而生了。
二.实验原理 1.FFT 的原理和参数生成公式: 我们先设序列长度为N=2^L,L 为整数。将N=2^L 的序列 x(n)(n=0,1,……,N-1),按N 的奇偶分成两组,也就是说我们将一个N 点的DFT 分解成两个N/2 点的DFT,他们又重新组合成一个如下式所表达的N 点DFT: 我们称这样的RFFT 优化算法是包装算法:首先2N 点实数的连续输入称为“进包”。其次N 点的FFT 被连续运行。最后作为结果产生的N 点的合成输出是“打开”成为最初的与DFT 相符合的2N 点输入。
三.FFT 的基本结构: 1.FFT 信号流图如下: 整个过程共有log2N 次,每次分组间隔为2^(L-1)----------------1=