第一篇:OFDM_matlab源程序总结
比较完整的OFDM仿真,供大家学习下载。、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、、%main_OFDM.m %一个相对完整的OFDM通信系统的仿真设计,包括编码,调制,IFFT,%上下变频,高斯信道建模,FFT,PAPR抑制,各种同步,解调和解码等模 %块,并统括系统性能的仿真验证了系统设计的可靠性。
clear all close all clc
%++++++++++++++++++++++++++全局变量++++++++++++++++++++++++++++++ % seq_num
表示当前帧是第几帧 % count_dds_up
上变频处的控制字的累加
% count_dds_down
下变频处的控制字的累加(整整)% count_dds_down_tmp
下变频处的控制字的累加(小数)% dingshi
定时同步的定位
% m_syn
记录定时同步中的自相关平台 global seq_num global count_dds_up global count_dds_down global count_dds_down_tmp global dingshi
global m_syn %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% SNR_Pre
设定用于仿真的信噪比的初值 % interval_SNR
设定用于仿真的信噪比间隔
% frame_num
每一个信噪比下仿真的数据帧数 % err_int_final
用于计算每帧出现的误比特数
% fwc_down
设定的接收机初始载波频率控制字
% fre_offset
设定接收机初始载波频率偏移调整量(单位为Hz)% k0
每次进入卷积编码器的信息比特数 % G
卷积编码的生成矩阵 SNR_Pre=-5;interval_SNR=1;
for SNR_System=SNR_Pre:interval_SNR:5
frame_num=152;dingshi=250;err_int_final=0;fwc_down=16.050;fre_offset=0;k0=1;G=[1 0 1 1 0 1 1;1 1 1 1 0 0 1 ];
disp('--------------start-------------------');
for seq_num=1:frame_num, %frame_num 帧数
%+++++++++++++++++++++++以下为输入数据部分+++++++++++++++++++++++ datain=randint(1,90);
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++以下为信道卷积编码部分+++++++++++++++++++++ encodeDATA=cnv_encd(G,k0,datain);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++++信道交织编码+++++++++++++++++++++++++ interlacedata=interlacecode(encodeDATA,8,24);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++以下为QPSK调制部分+++++++++++++++++++++ QPSKdata=qpsk(interlacedata);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++++生成训练序列+++++++++++++++++++++++++ if seq_num<3 trainsp_temp=seq_train();
end %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++++插入导频++++++++++++++++++++++++++++ PILOT=(1+j);m_QPSKdata=QPSKdata;data2fft_temp=[m_QPSKdata(1:8),PILOT,m_QPSKdata(9:16),PILOT,m_QPSKdata(17:24),PILOT,m_QPSKdata(25:32),PILOT,m_QPSKdata(33:40),PILOT,m_QPSKdata(41:48),m_QPSKdata(49:56),PILOT,m_QPSKdata(57:64),PILOT,m_QPSKdata(65:72),PILOT,m_QPSKdata(73:80),PILOT,m_QPSKdata(81:88),PILOT,m_QPSKdata(89:end)];%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
trainsp_temp2=[trainsp_temp,zeros(1,128)];trainsp=[trainsp_temp2(65:256),trainsp_temp2(1:64)];
%++++++++++++++++++++++++++降PAPR矩阵变换++++++++++++++++++++++++ matix_data=nyquistimp_PS();matrix_mult=data2fft_temp*matix_data;%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
data2fft2=[matrix_mult(65:128),zeros(1,128),matrix_mult(1:64)];
%++++++++++++++++++++++++++++ifft运算+++++++++++++++++++++++++++ if seq_num==1
ifftin=trainsp;elseif seq_num==2
ifftin=trainsp;else
ifftin=data2fft2;end IFFTdata=fft_my(conj(ifftin)/256);IFFTdata=conj(IFFTdata);% figure % plot(real(IFFTdata))% xlabel('realIFFTdata')% figure % plot(imag(IFFTdata))% xlabel('imagIFFTdata')
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++以下为插入循环前后缀,2倍升采样+++++++++++++++ data2fir=add_CYC_upsample(IFFTdata,2);% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% +++++++++++++++++++++++++fir低通滤波++++++++++++++++++++++++++ guiyi_a=[0.0017216 0.010162 0.025512 0.028801-0.0059219-0.060115-0.0496 0.091431 0.29636 0.3956 0.29636 0.091431-0.0496-0.060115-0.0059219 0.028801 0.025512 0.010162 0.0017216 ];%抽样截止频率为128kHZ,通带截止频率为20kHZ,阻带截止频率为40kHZ,带内纹波动小于1dB,带外衰减100dB txFIRdatai=filter(guiyi_a,1,real(data2fir));txFIRdataq=filter(guiyi_a,1,imag(data2fir));% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++发射机cic滤波++++++++++++++++++++++++++ CICidatai=cic_inter(txFIRdatai,20);CICidataq=cic_inter(txFIRdataq,20);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++++++上变频++++++++++++++++++++++++++++ fwc_up=16;
%控制字可以选择 DUCdata=up_convert_ofdm(fwc_up,CICidatai,CICidataq);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++高斯白噪声信道++++++++++++++++++++++++ [DUCdata,datamax]=guiyi_DUCdata(DUCdata);awgn_data=awgn(DUCdata,SNR_System);
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%*****************************接受机***************************** %+++++++++++++++++++++++++++++下变频+++++++++++++++++++++++++++++ DUCdata_tmp=awgn_data;fwc_down=fwc_down+(fre_offset*128/2560000);r_fre_offset=2560000*((fwc_down-fwc_up)/128);[DDCdatai,DDCdataq]=down_convert_ofdm(fwc_down,DUCdata_tmp);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++接收机cic滤波+++++++++++++++++++++++++ CICddatai=cic_deci(DDCdatai,40,40);CICddataq=cic_deci(DDCdataq,40,40);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++++fir低通滤波++++++++++++++++++++++++ guiyi_b=[ 0.019527-0.03934 0.049055-0.018102-0.1003 0.5944 0.5944-0.1003-0.018102 0.049055-0.03934 0.019527];%抽样截止频率为64kHZ,通带截止频率为20kHZ,阻带截止频率为30kHZ,带内纹波动小于1dB,带外衰减60dB rxFIRdatai=filter(guiyi_b,1,CICddatai);rxFIRdataq=filter(guiyi_b,1,CICddataq);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++++++++量化++++++++++++++++++++++++++++ q_rxFIRdatai=sign(rxFIRdatai);q_rxFIRdataq=sign(rxFIRdataq);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++++定时同步检测++++++++++++++++++++++++ if seq_num<3 time_syn(q_rxFIRdatai,q_rxFIRdataq);end %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++++++频率同步+++++++++++++++++++++++++++ fre_offset=fre_syn(rxFIRdatai,rxFIRdataq);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++++++fft运算+++++++++++++++++++++++++++ if seq_num>2
seq_num-2
fftw=32+dingshi;
rxFIRdata_syn=rxFIRdatai(fftw:fftw+255)+j*rxFIRdataq(fftw:fftw+255);
FFTdata=fft_my(rxFIRdata_syn);
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++降PAPR逆矩阵变换++++++++++++++++++++++ fftdata_reg=[FFTdata(193:256),FFTdata(1:64)];dematrix_data=fftdata_reg*pinv(matix_data);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++++++++++相位补偿+++++++++++++++++++++++++++ rx_qpsk_din_th=phase_comp(dematrix_data);%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++++QPSK解调部分++++++++++++++++++++++++ % figure % plot(rx_qpsk_din_th,'.')% xlabel('星座图')datatemp4=deqpsk(rx_qpsk_din_th);datatemp4=sign(datatemp4);for m=1:192
if datatemp4(m)==-1
datatemp4(m)=1;
elseif datatemp4(m)==1
datatemp4(m)=0;
end end %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++信道解交织++++++++++++++++++++++++++++ interdout=interlacedecode(datatemp4,8,24);%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%++++++++++++++++++++以下为viterbi译码部分++++++++++++++++++++++ decodeDATA=viterbi(G,k0,interdout);%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%+++++++++++++++++++++++++误比特统计++++++++++++++++++++++++++++ err_final=sum(abs(decodeDATA-datain))err_int_final=err_int_final+err_final end end disp('-----------------------');SNR_System err_rate_final((SNR_System-SNR_Pre)./interval_SNR+1)=err_int_final/(90*(frame_num-2))disp('-----------------------');end
disp('--------------end-------------------');
SNR_System=SNR_Pre:interval_SNR:5;figure semilogy(SNR_System,err_rate_final,'b-*');xlabel('信噪比/dB')ylabel('误码率')axis([-5,5,0,1])grid on %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%************************beginning of file***************************** %cnv_encd.m %卷积码编码程序
function output=cnv_encd(G,k0,input)% cnv_encd(G,k0,input),k0 是每一时钟周期输入编码器的 bit 数,% G 是决定输入序列的生成矩阵,它有 n0 行 L*k0 列 n0 是输出 bit 数,% 参数 n0 和 L 由生成矩阵 G 导出,L 是约束长度。L 之所以叫约束长度 % 是因为编码器在每一时刻里输出序列不但与当前输入序列有关,% 而且还与编码器的状态有关,这个状态是由编码器的前(L-1)k0。% 个输入决定的,通常卷积码表示为(n0,k0,m),m=(L-1)*k0 是编码 % 器中的编码存贮个数,也就是分为 L-1 段,每段 k0 个
% 有些人将 m=L*k0 定义为约束长度,有的人定义为 m=(L-1)*k0 % 查看是否需要补 0,输入 input 必须是 k0 的整数部
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++ % G
决定输入序列的生成矩阵
% k0
每一时钟周期输入编码器的 bit 数 % input
输入数据 % output
输入数据
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if rem(length(input),k0)>0 input=[input,zeros(size(1:k0-rem(length(input),k0)))];end n=length(input)/k0;% 检查生成矩阵 G 的维数是否和 k0 一致 if rem(size(G,2),k0)>0 error('Error,G is not of the right size.')end % 得到约束长度 L 和输出比特数 n0 L=size(G,2)/k0;n0=size(G,1);% 在信息前后加 0,使存贮器归 0,加 0 个数为(L-1)*k0 个 u=[zeros(size(1:(L-1)*k0)),input,zeros(size(1:(L-1)*k0))];% 得到 uu 矩阵,它的各列是编码器各个存贮器在各时钟周期的内容 u1=u(L*k0:-1:1);%将加 0 后的输入序列按每组 L*k0 个分组,分组是按 k0 比特增加 %从 1 到 L*k0 比特为第一组,从 1+k0 到 L*k0+k0 为第二组。。,%并将分组按倒序排列。for i=1:n+L-2 u1=[u1,u((i+L)*k0:-1:i*k0+1)];end uu=reshape(u1,L*k0,n+L-1);% 得到输出,输出由生成矩阵 G*uu 得到 output=reshape(rem(G*uu,2),1,n0*(L+n-1));% ************************end of file***********************************
%************************beginning of file***************************** %interlacecode.m
function dout=interlacecode(din,m,n)%实现信道的交织编码
%din为输入交织编码器的数据,m,n分别为交织器的行列值
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++ % din
输入数据
% m
交织器的行值 % n
交织器的列值 % dout
输出数据
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
for j=1:m
temp(j,:)=din(j*n-(n-1):j*n);end
dout_temp=reshape(temp,1,length(din));dout=dout_temp(1:end);%************************end of file**********************************
%************************beginning of file***************************** %qpsk.m %QPSK调制映射
function dout=qpsk(din)
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++ % din
输入数据 % dout
输出数据
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
din2=1-2*din;din_temp=reshape(din2,2,length(din)/2);for i=1:length(din)/2,dout(i)=din_temp(1,i)+j*din_temp(2,i);end % ************************end of file***********************************
%************************beginning of file***************************** %seq_train.m %生成用于同步的训练符号 function dout=seq_train()%第一帧产生短训练序列,第二帧产生长训练序列 %每个短训练符号由16个子载波组成,短训练序列 %是由伪随机序列经过数字调制后插0后,再经过 %IFFT之后得到的。具体过程如下:首先采用抽头 %系数为[1 0 0 1 ]的4级移位寄存器产生长度为
%15的伪随机序列之后末尾补0,经过QPSK调制之 %后的伪随机序列只在16的整数倍位置上出现,其 %余的位置补0,产生长度为128的序列,此序列再 %补128个0经过数据搬移后做256点的IFFT变换就 %得到16个以16为循环的训练序列,经过加循环前 %后缀就会产生20个相同的短训练序列。长训练序 %列的产生同短训练序列。
global seq_num
if seq_num==1
fbconnection=[1 0 0 1];
QPSKdata_pn=[m_sequence(fbconnection),0];
QPSKdata_pn=qpsk(QPSKdata_pn);elseif seq_num==2
fbconnection=[1 0 0 0 0 0 1];
QPSKdata_pn=[m_sequence(fbconnection),0];
QPSKdata_pn=qpsk(QPSKdata_pn);end
countmod=0;for k=1:128
if seq_num==1
if mod(k-1,16)==0
%生成16位循环的短训练符号
countmod=countmod+1;
trainsp_temp(k)=QPSKdata_pn(countmod);
else
trainsp_temp(k)=0;
end
elseif seq_num==2
if mod(k-1,2)==0
countmod=countmod+1;
trainsp_temp(k)=QPSKdata_pn(countmod);
else
trainsp_temp(k)=0;
end
end end
dout=trainsp_temp;% ************************end of file***********************************
%************************beginning of file***************************** %m_sequence.m %用线性移位寄存器产生m序列
function [mseq]= m_sequence(fbconnection);
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++ % fbconnection
线性移位寄存器的系数 % mseq
生成的m序列
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
n = length(fbconnection);N = 2^n-1;register = [zeros(1,n1.0000i
1.00001.0000i 1.00001.0000i
1.00001.0000i-1.00001.0000i
1.00001.0000i 1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.00001.0000i-1.00001.0000i
1.00001.0000i
1.0000 + 1.0000i
1.0000 + 1.0000i
1.0000 + 1.0000i
1.0000 + 1.0000i 1.0000 + 1.0000i-1.0000 + 1.0000i-1.00001.0000i
1.00001.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i
1.0000 + 1.0000i
1.0000 + 1.0000i-1.00001.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i
1.0000 + 1.0000i-1.0000 + 1.0000i-1.00001.0000i
1.0000 + 1.0000i
1.0000 + 1.0000i
1.0000 + 1.0000i-1.0000 + 1.0000i
1.0000 + 1.0000i
1.0000 + 1.0000i
1.00001.0000i-1.00001.0000i-1.00001.0000i-1.00001.0000i
1.00001.0000i
1.00001.0000i-1.00001.0000i
1.00001.0000i
1.00001.0000i
1.00001.0000i
1.00001.0000i-1.00001.0000i-1.00001.0000i-1.0000 + 1.0000i
1.0000 + 1.0000i
1.0000 + 1.0000i-1.0000 + 1.0000i-1.00001.0000i-1.00001.0000i-1.00001.0000i 1.0000-1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i-1.0000 + 1.0000i ];
for nn=1:length(datai)-128 %计算相关值(输入信号与本地信号互相关)
for m=1:100
t1_syn(m)=(datai(nn+m-1)+j*dataq(nn+m-1))*conj(local_seq(m));
end
lolol(nn)=sum(t1_syn);
t_syn(nn)=abs(sum(t1_syn));%输入信号与本地信号互相关判决函数
end
for ni=1:length(t_syn)
if t_syn(ni)>60
dingshi=find(t_syn(ni:ni+6)==max(t_syn(ni:ni+6)))+ni-1;
break
end end % figure % plot(t_syn)
% xlabel('采样点索引号')% ylabel('互相关判决函数')end % ************************end of file***********************************
%************************beginning of file***************************** %fre_syn.m
function dout=fre_syn(datai,dataq)%实现系统的频率同步
%频偏的估计范围由帧长度和循环间隔长度来决定。在 %本系统中,短训练序列的频偏估计范围为(2000Hz);长训 %练序列频偏估计范围为250Hz,前一个有较大的纠偏范围,估 %计值得到的方差较大;后一个方法的纠偏范围较小,估计值得 %到的方差较小。频率跟踪可以用循环前后缀的周期重复性来完 %成,其频偏估计范围为125Hz。
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++ % datai
输入数据的实部 % dataq
输入数据的虚部 % dout
输出数据
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
global seq_num
global dingshi
global m_syn
if seq_num==1
if m_syn(50)>10
for m=1:128
fre_back_tmp(m)=(datai(40+m-1)+j*dataq(40+m-1))*conj(datai(40+m-1+16)+j*dataq(40+m-1+16));
end
fre_back_sum=sum(fre_back_tmp);
fre_offset_tmp=-320*angle(fre_back_sum)/(2*pi*16);
fre_offset=fre_offset_tmp/0.005;
end elseif seq_num==2
for m=1:128
fre_back_tmp(m)=(datai(dingshi+m-1)+j*dataq(dingshi+m-1))*conj(datai(dingshi+m-1+128)+j*dataq(dingshi+m-1+128));
end
fre_back_sum=sum(fre_back_tmp);
fre_offset_tmp=-320*angle(fre_back_sum)/(2*pi*128);
fre_offset=fre_offset_tmp/0.005;
elseif seq_num>2
for m=1:48
fre_back_tmp(m)=(datai(dingshi+m)+j*dataq(dingshi+m))*conj(datai(dingshi+m+256)+j*dataq(dingshi+m+256));
end
fre_back_sum=sum(fre_back_tmp);
fre_offset_tmp=-320*angle(fre_back_sum)/(2*pi*256);
fre_offset= fre_offset_tmp/0.005;end dout=fre_offset;% ************************end of file**********************************
%************************beginning of file***************************** %phase_comp.m %实现OFDM符号的相位补偿
function
dout=phase_comp(din)
%使用导频的相位信息来调整有效数据的偏移相位,%以实现正确的QPSK解调
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++ % din
输入数据 % dout
输出数据
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
qpsk_din=[din(1:8),din(10:17),din(19:26),din(28:35),din(37:44),din(46:53),din(54:61),din(63:70),din(72:79),din(81:88),din(90:97),din(99:106)];pilot=[din(9),din(18),din(27),din(36),din(45),din(62),din(71),din(80),din(89),din(98)];ang_offset=angle(sum((1-j)*pilot));qpsk_din=qpsk_din*exp(-j*ang_offset);dout=qpsk_din;% ************************end of file***********************************
%************************beginning of file***************************** %deqpsk.m
function dout=deqpsk(din)%实现QPSK解调
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++ % din
输入数据 % dout
输出数据
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
dout_temp(1,:)=real(din);dout_temp(2,:)=imag(din);
dout=reshape(dout_temp,1,length(din)*2);% ************************end of file**********************************
%************************beginning of file***************************** %interlacedecode.m
function dout=interlacedecode(din,m,n)%实现信道的交织解码
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++ % din
输入数据
% m
交织器的行值 % n
交织器的列值 % dout
输出数据
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
temp=reshape(din,m,n);for j=1:m
dout_temp(j*n-(n-1):j*n)=temp(j,:);end
dout=dout_temp(1:end);% ************************end of file**********************************
%************************beginning of file***************************** %viterbi.m %viterbi 解码程序
% function [decoder_output,survivor_state,cumulated_metric]=viterbi(G,k,channel_output)function decoder_output=viterbi(G,k,channel_output)% [decoder_output,survivor_state,cumulated_metric]=viterbi(G,k,channel_output)% 其中 G 是一个 n 行 L*k 列矩阵,它的每一行决定了从移位寄存器到输入码字的连接方式.% survivor_state 是一个矩阵,它显示了通过网格的最优路径,这个矩阵通过一个单独 % 的函数 metric(x,y)给出。
n=size(G,1);%检验 G 的维数 if rem(size(G,2),k)~=0 error('Size of G and k do not agree')end if rem(size(channel_output,2),n)~=0 error('channle output not of the right size')end L=size(G,2)/k;number_of_states=2^((L-1)*k);%产生状态转移矩阵,输出矩阵和输入矩阵 for j=0:number_of_states-1 for t=0:2^k-1 [next_state,memory_contents]=nxt_stat(j,t,L,k);input(j+1,next_state+1)=t;branch_output=rem(memory_contents*G',2);nextstate(j+1,t+1)=next_state;output(j+1,t+1)=bin2deci(branch_output);end end input;state_metric=zeros(number_of_states,2);depth_of_trellis=length(channel_output)/n;channel_output_matrix=reshape(channel_output,n,depth_of_trellis);survivor_state=zeros(number_of_states,depth_of_trellis+1);[row_survivor col_survivor]=size(survivor_state);%开始非尾信道输出的解码
%i 为段,j 为每一阶段的状态,t 为输入 for i=1:depth_of_trellis-L+1 flag=zeros(1,number_of_states);if i<=L step=2^((L-i)*k);else step=1;end for j=0:step:number_of_states-1 for t=0:2^k-1 branch_metric=0;binary_output=deci2bin(output(j+1,t+1),n);for tt=1:n branch_metric=branch_metric+metric(channel_output_matrix(tt,i),binary_output(tt));end if((state_metric(nextstate(j+1,t+1)+1,2)>state_metric(j+1,1)+branch_metric)|flag(nextstate(j+1,t+1)+1)==0)state_metric(nextstate(j+1,t+1)+1,2)=state_metric(j+1,1)+branch_metric;survivor_state(nextstate(j+1,t+1)+1,i+1)=j;flag(nextstate(j+1,t+1)+1)=1;end end end state_metric=state_metric(:,2:-1:1);end %开始尾信道输出的解码
for i=depth_of_trellis-L+2:depth_of_trellis flag=zeros(1,number_of_states);last_stop=number_of_states/(2^((i-depth_of_trellis+L-2)*k));for j=0:last_stop-1 branch_metric=0;binary_output=deci2bin(output(j+1,1),n);for tt=1:n branch_metric=branch_metric+metric(channel_output_matrix(tt,i),binary_output(tt));end if((state_metric(nextstate(j+1,1)+1,2)>state_metric(j+1,1)+branch_metric)|flag(nextstate(j+1,1)+1)==0)state_metric(nextstate(j+1,1)+1,2)=state_metric(j+1,1)+branch_metric;survivor_state(nextstate(j+1,1)+1,i+1)=j;flag(nextstate(j+1,1)+1)=1;end end state_metric=state_metric(:,2:-1:1);end %从最优路径产生解码输出 %由段得到状态序列,再由状序列从 input 矩阵中得到该段的输出 state_sequence=zeros(1,depth_of_trellis+1);size(state_sequence);state_sequence(1,depth_of_trellis)=survivor_state(1,depth_of_trellis+1);for i=1:depth_of_trellis state_sequence(1,depth_of_trellis-i+1)=survivor_state((state_sequence(1,depth_of_trellis+2-i)+1),depth_of_trellis-i+2);end state_sequence;decoder_output_matrix=zeros(k,depth_of_trellis-L+1);for i=1:depth_of_trellis-L+1 dec_output_deci=input(state_sequence(1,i)+1,state_sequence(1,i+1)+1);dec_output_bin=deci2bin(dec_output_deci,k);decoder_output_matrix(:,i)=dec_output_bin(k:-1:1)';end decoder_output=reshape(decoder_output_matrix,1,k*(depth_of_trellis-L+1));cumulated_metric=state_metric(1,1);% ************************end of file***********************************
%************************beginning of file***************************** %nxt_stat.m %viterbi 解码子程序
function [next_state,memory_contents]=nxt_stat(current_state,input,L,k)
binary_state=deci2bin(current_state,k*(L-1));binary_input=deci2bin(input,k);next_state_binary=[binary_input,binary_state(1:(L-2)*k)];next_state=bin2deci(next_state_binary);memory_contents=[binary_input,binary_state];%************************end of file***********************************
%************************beginning of file***************************** %deci2bin.m function y=deci2bin(x,t)% 十进制数 x 转化为二进制数,二进制数至少表示为 t 位
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++ % x
输入数据 % t
二进制数位数 % y
输出数据
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
y=zeros(1,t);i=1;while x>=0&i<=t y(i)=rem(x,2);x=(x-y(i))/2;i=i+1;end y=y(t:-1:1);% ************************end of file**********************************
%************************beginning of file***************************** %bin2deci.m function y=bin2deci(x)% 将二进制数转化为十进制数
%+++++++++++++++++++++++variables++++++++++++++++++++++++++++ % x
输入数据 % y
输出数据
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
t=length(x);y=(t-1:-1:0);y=2.^y;y=x*y';% ************************end of file*********************************
%************************beginning of file***************************** %metric.m %viterbi 解码子程序
function distance=metric(x,y)if x==y distance=0;else distance=1;end % switch(y)% case 0 % if x==0 % distance=0.0458;% end % if x==1 % distance=2;% end % if x==2;% distance==1.0458;% end % case 1 % if x==0 % distance==2;% end % if x==1 % distance=0.0458;% end % if x==2 % distance=1.0458;% end % otherwise, % break;% end % ************************end of file**********************************
%完结
第二篇:病房管理源程序
附件:病房管理源程序
#include
typedef struct bed{
int
number;
/*床位号*/ int
live;
/*是否有人*/
unsigned char name[max_name];
/*人名最长长度*/ struct bed
*next;
/*下一个床铺*/ }bed;
typedef struct sickroom { int
number;
/*病房号*/ int
bednum;
/*房内床位数*/ int
livenum;
/*已住人数*/ struct sickroom *next;
/*指向下一个病房*/
struct bed *firstbed;
/*指向第一个床位*/ }sickroom;
int roomnum(sickroom *q)/*返回医院中的房间的个数*/ {
int i;
for(i=0;q!=NULL;i++)q=q->next;return i;}
int Initialization(sickroom *H)/*初始化操作*/ {
int i,n,k;sickroom *q,*sq;bed *p,*bp;printf(“----------------------请输入病房个数:”);scanf(“%d”,&n);getchar();
int i,x,n,m,y;sickroom *q;bed *p;q=H->next;n=roomnum(q);
printf(“请输入您要入住的病房号(目前共有%d个病房):”,n);scanf(“%d”,&x);getchar();printf(“n”);while(x>n||x<1){
printf(“您输入的病房号错误,请重新输入:”);
scanf(“%d”,&x);
getchar();
printf(“n”);} for(i=1;i
/*q指向选择的病房*/ m=q->bednum-q->livenum;if(m==0){
printf(“-------------------------此病房已满。。(回车继续)-----------n”);
getchar();
printf(“n”);
return 0;} p=q->firstbed;printf(“空床铺有:”);for(i=1;i<=q->bednum;i++){
if(p->live==0)printf(“ %d号床铺”,i);
p=p->next;} printf(“n输入您要选择的床铺号:”);scanf(“%d”,&y);getchar();p=q->firstbed;for(i=1;i
/*指向所选床铺*/ if(p->live==1){
printf(“-------------------------输入错误!(回车继续)----------------------------n”);
getchar();
return 0;} p->live=1;
{
printf(“-------------------------查无此床(回车继续)-----------n”);
getchar();
return 0;
} } if(q->live==0){
printf(“-------------------------此床铺为空!(回车继续)-------------n”);
getchar();
return 0;} q->live=0;p->livenum--;H->livenum--;printf(“-------------------------”);printf(“%s”,q->name);printf(“已成功出院!n”);printf(“------------------------(回车继续)----------------n”);getchar();for(i=1;i<=21;i++)q->name[i-1]=0;return 1;}
int searchsickroom(sickroom *H)/*查找出医院内病房空床位*/ {
int i,m,z=0;
sickroom *q;
q=H->next;
for(i=1;q;i++)
{
m=q->bednum-q->livenum;
if(m>0)printf(“
%d号病房有%d个床位n”,i,m);
else printf(“
%d号病房已满n”,i);
q=q->next;
z+=m;
}
printf(“
院中还剩%d个床位(回车继续)。。n”,z);
getchar();
return 1;}
int n=0,m=0;
sickroom *q;
bed *p;
char a[max_name];
printf(“n请输入您要查找的人名(二十个字符以内,以空格结束):”);
scanf(“%s”,a);
getchar();
for(q=H->next;q;q=q->next)
{
for(p=q->firstbed;p;p=p->next)
{
if(p->live==1)
if(strcmp(a,p->name)==0)
{
n=q->number;
m=p->number;
}
}
}
if(m==0&&n==0)
{
printf(“-------------------------查无此人(回车继续)------------n”);
getchar();
}
else
{
printf(“
此人目前住在%d号房%d号床。。n”,n,m);
printf(“---------------------(回车继续)-----------------n”);
getchar();
}
return 1;}
int nowhospital(sickroom *H)/*显示目前医院的住院情况*/ {
int i;
sickroom *q;
bed *p;
i=0;
for(q=H->next;q;q=q->next)
{
printf(“n病房%d —>t”,q->number);
for(p=q->firstbed;p;p=p->next)
if(i==2)searchbed(H);
else
if(i==3)
searchpeople(H);
else
if(i==4)
nowhospital(H);
else
if(i==0);
else
{
printf(“输入错误!(回车继续)n”);
getchar();
}
}while(i!=0);
return ok;}
void main()/*主程序*/ {
int i;
sickroom H;
system(“CLS”);
printf(“********12**************14************12************14********12*********12 n”);
printf(“
_______________________________________________________
n”);
printf(“
n”);
printf(“
欢迎进入医院管理程序!!
n”);
printf(“
n”);
printf(“
请先对医院进行初始化。。(回车键继续)
n”);
printf(“
1_____________________________________________________ 2
n”);
getchar();
Initialization(&H);
printf(“n”);
do
{
system(“CLS”);
printf(“********12**************14************12************14********12******12********n”);
printf(“
_____________________________________
n”);
printf(“
n”);
printf(“
选择菜单:
n”);
printf(“
[1] 入院操作
n”);
printf(“
[2] 出院操作
n”);
printf(“
[3] 查询医院信息
n”);
printf(“
[0] 退出医院
n”);
第三篇:c语言源程序
基于单片机msp430和温度传感器ds18b20的水温度控制系统的c语言源程序(不是测量,要有加热跟制冷)
我这是用STC做的,应该很容易移植到MPS430上的给你参考一下。#include
sbit scl=P1^3;sbit sda=P1^4;
sbit key1=P1^6;sbit key2=P1^7;sbit key3=P2^0;sbit key4=P2^1;
sbit lcrs=P3^7;//数据/命令 sbit lcwr=P3^5;//读/写 sbit lcden=P3^4;//使能
sbit DS=P2^2;
/*sbit lcrs=P3^4;//数据/命令 sbit lcwr=P3^7;//读/写 sbit lcden=P3^5;//使能 */ sbit jrk=P2^2;sbit cyk=P2^3;sbit xhk=P2^4;bit flag=0,rsg=0,not=0,he=0,in=0;int acon=0,bcon=0,dcon=0,econ=0, temp=0,y=0,j=0,l=0,cfj=0,ec=0,dc=0,at;uchar code table[]={48,49,50,51,52,53,54,55,56,57};uchar code ta1[]={“Temperature UP”};uchar code ta2[]={“Temperature DN”};uchar code ta3[]={“Inflator Cycle”};uchar code ta4[]={“Inflator Time ”};uchar code ta5[]={“ Heating UP ”};uchar code ta6[]={“ Inflator ”};uchar code table7[]={“Temperature”};uchar table1[]={0,0,0,'.',0};uchar table3[]={“AptitudeAquarium”};uchar table4[]={0,0,0,0,0};uchar n,c=0;void delay(uchar);void wen_kong();void xh();void rso();void weno();
void Init_Com(void){ TMOD = 0x11;PCON = 0x00;TH1=0x61;TL1=0x99;EA=1;ET1=1;TR1=1;} void delay(uchar count)//delay { uint i;while(count){ i=200;while(i>0)i--;count--;} } ////初始化18B20///////// bit init18b20(void){ uint i;bit no;DS=0;i=103;while(i>0)i--;DS=1;i=4;while(i>0)i--;no=DS;if(no==0){ DS=1;i=100;while(i>0)i--;no=DS;if(no==1)not=0;else not=1;} else not=1;return(not);}
bit tmpread()bit(void)//读一位
{ uint i;bit dat;DS=0;i++;DS=1;i++;i++;dat=DS;i=8;while(i>0)i--;return(dat);}
uchar tmpread()(void)//读一个字节 { uchar i,j,dat;dat=0;for(i=1;i<=8;i++){ j=tmpread()bit();dat=(j<<7)|(dat>>1);//读出的数据最低位在最前面,这样刚好一个字节在DAT里
} return(dat);}
void tmpwritebyte(uchar dat)//写一个字节到 ds18b20 { uint i;uchar j;bit testb;for(j=1;j<=8;j++){ testb=dat&0x01;dat=dat>>1;if(testb)//write 1 { DS=0;i++;i++;DS=1;i=8;while(i>0)i--;} else { DS=0;//write 0 i=8;while(i>0)i--;DS=1;i++;i++;} } }
int tmp()//DS18B20温度读取 { float tt;int a,b;if(init18b20()==0){ WDT_CONTR=0x36;/////喂狗 EA=0;delay(1);
tmpwritebyte(0xcc);// 跳过读ROM操作 tmpwritebyte(0x44);// 启动温度转换 delay(10);init18b20();delay(1);tmpwritebyte(0xcc);tmpwritebyte(0xbe);a=tmpread();b=tmpread();temp=b;temp<<=8;//将高字节温度数据与低字节温度数据整合 temp=temp|a;c=b>>4;
tt=temp*0.0625;
temp=tt*10+0.5;//放大10倍输出并四舍五入 EA=1;return temp;} else not=1;}
//////1062///////// void ydelay(uint x){ uint a,b;for(a=x;a>0;a--)for(b=10;b>0;b--);} void write_com(uchar com){ P0=com;lcwr=0;lcrs=0;lcden=0;ydelay(10);lcden=1;ydelay(10);lcden=0;lcwr=1;}
void write_date(uchar date)//写数据 {
P0=date;lcwr=0;lcrs=1;lcden=0;ydelay(10);lcden=1;ydelay(10);lcden=0;lcwr=1;}
void init1602()//初始化 { write_com(0x38);//设置显示模式 ydelay(20);write_com(0x0c);//开显示 ydelay(20);write_com(0x06);//指针和光标自动加一 ydelay(20);write_com(0x01);//清屏指令 ydelay(20);}
///////显示程序////// void display(int num){ uint i,A1,A2;WDT_CONTR=0x35;/////喂狗 if(c!=0)num=~num+1;A1=num/1000;A2=num%1000/100;if(not==0){
if(c!=0){ c=0;table1[0]='-';} else if(A1==0)table1[0]=' ';else
table1[0]=table[A1];if(A1==0)if(A2==0)table1[1]=' ';else table1[1]=table[A2];
table1[2]=table[num%1000%100/10];table1[4]=table[num%1000%100%10];} else { table1[0]='?';table1[1]='?';table1[2]='?';table1[4]='?';}
write_com(0x80);for(i=0;i<11;i++){write_date(table7[i]);delay(2);} write_com(0x8b);for(i=0;i<5;i++){write_date(table1[i]);delay(2);} write_com(0xc0);for(i=0;i<16;i++){ if(he==1)write_date(ta5[i]);else if(in==1)write_date(ta6[i]);else write_date(table3[i]);} c=0;WDT_CONTR=0x35;/////喂狗 } ////显示2//////////////////// display2(uchar bh,int dat){ uchar a,A,B;WDT_CONTR=0x35;/////喂狗 //write_com(0x01);//清屏指令 y=dat;y=y&0x8000;if(y!=0)dat=~dat+1;A=dat/1000;B=dat%1000/100;if((bh!=4)&&(bh!=5)){ if(A!=0)table4[0]=table[dat/1000];else if((c!=0)||(y!=0)){ c=0;y=0;table4[0]='-';} else table4[0]=' ';if(B!=0)table4[1]=table[B];else table4[1]=' ';table4[2]=table[dat%1000%100/10];table4[3]='.';table4[4]=table[dat%1000%100%10];} else { table4[0]=' ';if((c!=0)||(y!=0)){ c=0;y=0;table4[1]='-';} else table4[1]=' ';table4[2]=' ';table4[3]=table[dat%1000%100/10];table4[4]=table[dat%1000%100%10];}
write_com(0xc4);delay(2);for(a=0;a<5;a++)write_date(table4[a]);delay(2);write_com(0x80);switch(bh){ case 1:for(a=0;a<14;a++)write_date(ta1[a]);break;case 2:for(a=0;a<14;a++)write_date(ta2[a]);break;case 3:for(a=0;a<14;a++)write_date(ta3[a]);break;case 4:for(a=0;a<14;a++)write_date(ta4[a]);break;default:break;} }
///////////x24c02////////////////// void delay24(){;;}
void init24c02()//初始化 { sda=1;delay24();scl=1;delay24();}
void start()//开始信号 { sda=1;delay24();scl=1;delay24();sda=0;delay24();}
void stop()//停止 { sda=0;delay24();scl=1;delay24();sda=1;delay24();}
void respons()//应答 { uchar i;scl=1;delay24();while((sda==1)&&(i<250))i++;scl=0;delay24();}
void write_byte(uchar date)// 写数据子函数 { uchar i,temp;temp=date;
for(i=0;i<8;i++){ temp=temp<<1;scl=0;delay24();sda=CY;delay24();scl=1;delay24();} scl=0;delay24();sda=1;delay24();}
uchar read_byte()// 读数据子函数 { uchar i,k;scl=0;delay24();sda=1;delay24();for(i=0;i<8;i++){ scl=1;delay24();k=(k<<1)|sda;scl=0;delay24();} return k;} ///////写数据函数/////////////////// void write_add(uchar address,uint date){ start();write_byte(0xa0);respons();write_byte(address);respons();write_byte(date/256);respons();write_byte(date%256);respons();stop();} uchar read_add(uchar address)//读数据函数 { uchar date;start();write_byte(0xa0);respons();write_byte(address);respons();start();write_byte(0xa1);respons();date=read_byte();stop();return date;}
void delay1ms(uchar ms){
uchar i;while(ms--){ for(i = 0;i< 250;i++){ _nop_();_nop_();_nop_();_nop_();} } }
int keyf(int *num,int up,int dn){ uint i;uchar z;for(i=0;i<600;i++){
display2(n,*num);if(key1==0){ delay1ms(30);if(key1==0){ i=0;n++;if(n>=9)n=0;while(!key1)display2(n,*num);break;} } if(key2==0){
delay1ms(10);if(key2==0){ i=0;if(*num>=up)*num=up;else if(n!=4)*num=*num+1;else if(*num<100)*num=*num+5;else
*num=*num+10;for(z=0;z<65;z++){ display2(n,*num);if(key2!=0)break;} while(!key2){ for(z=0;z<2;z++)display2(n,*num);if(*num>=up)*num=up;else if(n!=4)*num=*num+1;else if(*num<100)*num=*num+5;else
*num=*num+10;} } }
if(key3==0){ delay1ms(10);if(key3==0){ i=0;if(*num<=dn)*num=dn;else if(n!=4)*num=*num-1;else if(*num<100)*num=*num-5;else
*num=*num-10;for(z=0;z<65;z++){ display2(n,*num);if(key3!=0)break;} while(!key3){ for(z=0;z<2;z++)display2(n,*num);if(*num<=dn)*num=dn;else if(n!=4)*num=*num-1;else if(*num<100)*num=*num-5;else
*num=*num-10;} } } } return(*num);} void keyjc(){ uchar i=0;if(key1==0){ delay1ms(10);if(key1==0){ EA=0;
for(i=0;i<20;i++){
display(tmp());} if(key1==0){
write_com(0x01);//清屏指令
n++;if(n>=5)n=0;while(!key1){ switch(n){ case 1:display2(n,acon);break;case 0:break;} } if(n==1){ keyf(&acon,1250,-530);if((acon-bcon)<3)bcon=acon-3;} if(n==2){ keyf(&bcon,1240,-550);if((acon-bcon)<3)acon=bcon+3;} write_add(1,acon);//A delay1ms(15);write_add(3,bcon);//B n=0;write_com(0x01);//清屏指令 } EA=1;} } }
key(){ uint i;if(key4==0)delay1ms(50);if(key4==0){ write_com(0x01);//清屏指令
for(i=0;i<500;i++){ if(key4==0){ delay1ms(15);if(key4==0){ i=0;n++;if(n>=5)n=0;while(!key4){ switch(n){ case 1: display2(1,acon);break;case 2: display2(2,bcon);break;default: break;} } } } switch(n){ case 1: display2(1,acon);break;case 2: display2(2,bcon);break;default: break;} } n=0;} }
///////滤波//////// int filter(){ int tm,buf[6];uchar i,j;EA=0;for(i=0;i<6;i++){ buf[i]=tmp();delay1ms(20);WDT_CONTR=0x35;/////喂狗 }
for(j=0;j<5;j++)for(i=0;i<5-j;i++)if(buf[i]>buf[i+1]){ tm=buf[i];buf[i]=buf[i+1];buf[i+1]=tm;} tm=((buf[2]+buf[3])/2);EA=1;return(tm);}
void main(){ uchar b,c;Init_Com();init1602();init24c02();
b=read_add(1);delay1ms(15);c=read_add(2);delay1ms(15);acon=b*256+c;b=read_add(3);delay1ms(15);c=read_add(4);delay1ms(15);bcon=b*256+c;
AUXR=0x01;// 禁止ALE输出 WDT_CONTR=0x35;//启动看门狗 write_com(0x01);//清屏指令
while(1){ at=filter();display(at);keyjc();key();
wen_kong();weno();} }
//////温度控制//////////////
void wen_kong(){ if((flag==0)&&(not==0)){ at=filter();if(at<=bcon)
{ flag=1;jrk=0;xhk=0;he=1;} } }
void weno(){ if(flag){ at=filter();if(at>=acon){ flag=0;jrk=1;if(rsg)xhk=0;else xhk=1;he=0;} } if(not==1){ flag=0;jrk=1;if(rsg)xhk=0;else xhk=1;he=0;} }
第四篇:图像放大算法总结及MATLAB源程序
1,插值算法(3种):
(1)最邻近插值(近邻取样法):
最邻近插值的的思想很简单,就是把这个非整数坐标作一个四舍五入,取最近的整数点坐标处的点的颜色。可见,最邻近插值简单且直观,速度也最快,但得到的图像质量不高。
最邻近插值法的MATLAB源代码为:
A = imread('F:lena.jpg');%读取图像信息 imshow(A);%显示原图 title('原图128*128');
Row = size(A,1);Col = size(A,2);%图像行数和列数 nn=8;%放大倍数
m = round(nn*Row);%求出变换后的坐标的最大值 n = round(nn*Col);
B = zeros(m,n,3);%定义变换后的图像
for i = 1 : m
for j = 1 : n
x = round(i/nn);y = round(j/nn);%最小临近法对图像进行插值
if x==0 x = 1;end
if y==0 y = 1;end
if x>Row x = Row;end
if y>Col y = Col;end B(i,j,:)= A(x,y,:);
end end
B = uint8(B);%将矩阵转换成8位无符号整数 figure;imshow(B);
title('最邻近插值法放大8倍1024*1024');
运行程序后,原图如图1所示:
图1
用最邻近插值法放大4倍后的图如图2所示:
图2
(2)双线性内插值法:
在双线性内插值法中,对于一个目的像素,设置坐标通过反向变换得到的浮点坐标为(i+u,j+v),其中i、j均为非负整数,u、v为[0,1)区间的浮点数,则这个像素得值 f(i+u,j+v)可由原图像中坐标为(i,j)、(i+1,j)、(i,j+1)、(i+1,j+1)所对应的周围四个像素的值决定,即:
f(i+u,j+v)=(1-u)(1-v)f(i,j)+(1-u)vf(i,j+1)+ u(1-v)f(i+1,j)+ uvf(i+1,j+1)其中f(i,j)表示源图像(i,j)处的的像素值,以此类推。
这就是双线性内插值法。双线性内插值法计算量大,但缩放后图像质量高,不会出现像素值不连续的的情况。由于双线性插值具有低通滤波器的性质,使高频分量受损,所以可能会使图像轮廓在一定程度上变得模糊。
在MATLAB中,可用其自带的函数imresize()来实现双线性内插值算法。
双线性内插值算法的MATLAB源代码为:
A=imread('F:lena.jpg');imshow(A);
title('原图128*128');
C=imresize(A,8,'bilinear');%双线性插值 figure;imshow(C);
title('双线性内插值法放大8倍1024*1024');
程序运行后,原图如图3所示:
图3
双线性内插值法放大8倍后的图如图4所示:
图4
(3)双三次插值法:
双三次插值法能够在很大程度上克服以上两种算法的不足,该算法计算精度高,但计算量大,它考虑一个浮点坐标(i+u,j+v)周围的16个邻点。
目的像素值f(i+u,j+v)可由如下插值公式得到:f(i+u,j+v)= [A] * [B] * [C] 其中[A]=[ S(u + 1)S(u + 0)S(u2)];[C]=[ S(v + 1)S(v + 0)S(v2)];而[B]是周围16个邻点组成的4*4的矩阵;S(x)是对 Sin(x*π)/x 的逼近。
在MATLAB中,可用其自带的函数imresize()来实现双三次插值算法。MATLAB源代码为:
A=imread('F:lena.jpg');%读取原图像
D=imresize(A,8,'bicubic');%双三次插值放大8倍 figure;
T
imshow(D);title('三次卷积法放大8倍1024*1024');
MATLAB自带双三次插值法运行结果如图5所示:
图5
也可以自己编写双三次插值算法MATLAB代码如下:
clc,clear;
ff=imread('F:lena.jpg');%读取图像到ff
k=8;%设置放大倍数 [m,n,color]=size(ff);
f=zeros(m,n);%将彩色图像ff转换为黑白图像f for i=1:m
for j=1:n
f(i,j)=ff(i,j);
end end
a=f(1,:);c=f(m,:);%将待插值图像矩阵前后各扩展两行两列,共扩展四行四列 b=[f(1,1),f(1,1),f(:,1)',f(m,1),f(m,1)];d=[f(1,n),f(1,n),f(:,n)',f(m,n),f(m,n)];a1=[a;a;f;c;c];a1';
b1=[b;b;a1';d;d];f=b1';f1=double(f);
for i=1:k*m %利用双三次插值公式对新图象所有像素赋值 u=rem(i,k)/k;i1=floor(i/k)+2;A=[sw(1+u)sw(u)sw(1-u)sw(2-u)];
for j=1:k*n
v=rem(j,k)/k;j1=floor(j/k)+2;C=[sw(1+v);sw(v);sw(1-v);sw(2-v)];
B=[f1(i1-1,j1-1)f1(i1-1,j1)f1(i1-1,j1+1)f1(i1-1,j1+2)f1(i1,j1-1)f1(i1,j1)f1(i1,j1+1)f1(i1,j1+2)f1(i1,j1-1)f1(i1+1,j1)f1(i1+1,j1+1)f1(i1+1,j1+2)f1(i1+2,j1-1)f1(i1+2,j1)f1(i1+2,j1+1)f1(i1+2,j1+2)];g1(i,j)=(A*B*C);
end end
g=uint8(g1);%将矩阵转换成8位无符号整数 imshow(g);
title('自编双三次插值法放大8倍图像');
其中子函数sw代码如下: function A=sw(w1)w=abs(w1);if w<1&&w>=0 A=1-2*w^2+w^3;elseif w>=1&&w<2 A=4-8*w+5*w^2-w^3;else
A=0;end
与MATLAB自带函数相比,以上手工编写的MATLAB代码只能完成黑白图像输出,且运行时间远比MATLAB自带函数的运行时间长。手工编写双三次插值算法MATLAB代码的运行结果如图6所示:
图6
2,其他算法简介:
传统的图像放大方法有重复放大线性放大和高次多项式插值放大。重复放大最简单,但会产生明显的方块效应线性放大消除了方块效应,但会造成图像的模糊 高次多项式插值放大效果较好,但运算复杂。由于传统方法的固有缺陷,诞生了新一代图像放大方法,主要有小波放大、邻域交换内插和分形放大等。
下面简单介绍一下增强系数小波放大算法: 算法示意图如图7所示:
图7 通过二维离散小波变换,经分析高通滤波器和分析低通滤波器,可将一幅分辨率为p的二维图像分解为分辨率为p/2的离散逼近信号A1和水平、垂直、对角三个细节信号H1、V1、D1。这四个分量都只有原图像大小的1/4。之后又可以对A1进行同样的分解如图7所示。这个过程可以一直重复下去。通过二维离散小波反变换,用相应的综合高通滤波器和综合低通滤波器可将各分量重构为原图像。
对于一个图像,低频成分包含了基本特征,即原图像的近似,高频成分反应其细节。基于此,我们将原图像作为低频成分A1,其他3个细节部分置0,进行小波重构,便可得到放大4倍的图像。但是由于能量守恒,放大后的图像能量分散会显得较暗。可以将原图像灰度值矩阵乘2,再进行上述变换,便可解决这一问题。小波分解重构是一种全局运算,不会造成重复放大中的方块效应,同时较好地保持图像边缘的清晰。
第五篇:c语言源程序段
1.有三个整数a,b,c,由键盘输入,输出其中最大的数。#include
scanf(“%d%d%d”,&a,&b,&c);if(a>b&&a>c)
printf(“%dn”,a);else
if(b>a&&b>c)
printf(“%dn”,b);
else
printf(“%dn”,c);} 2.编程输入整数a和b,若a2b2大于100,则输出a2b2百位以上的数字,否则输出两数之和。
#include
printf(“%dn”,(d+e));else
printf(“%dn”,f);} 3.有一函数:
x(x1) y2x11(1x10)
3x11(x10)编写一程序,输入x,输出y值。#include
y=x;else
if(x<10&&x>=1)
y=2*x-11;
else
y=3*x-11;
printf(“%dn”,y);} 4.给出一百分制成绩,要求输出成绩等级’A’,’B’,’C’,’D’,’E’。90分以上为’A’,80-89分为’B’,70-79分为’C’,60-69分为’D’,60分以下为’E’ #include
printf(“A”);else
if(x<90&&x>=80)
printf(“B”);
else
if(x<80&&x>=70)
printf(“C”);
else
if(x<70&&x>=60)
printf(“D”);
else
printf(“E”);printf(“n”);} 5.提高题:给一个不多于5位的正整数,要求:①求出它是几位数;②分别打印出每一位数字;③按逆序打印出各位数字,例如原数是321,应输出123。#include
for(i=0;j>1;i++)
j=j/10;}
printf(“%dnn”,i);{
for(k=1;k<=i;k++){
b=a%10;
a=a/10;
printf(“%d”,b);} } }.求解一元二次方程a*x2+b*x+c=0 #include
int main(){ int a,b,c,m;double x1,x2,n;
//解为double类型
printf(“请输入ax2+bx+c=0中的a,b,c: n”);scanf(“%d %d %d”,&a,&b,&c);//输入参数
m=(b*b-4*a*c);if(m<0)
printf(“方程无解”);else
n=sqrt((double)m);
//对m进行强制类型转换为double,因为接为double
x1=(-b-m)/(2*(double)a);
x2=(-b+m)/(2*(double)a);
printf(“x1=%.2lf x2=%.2lfn”,x1,x2);return 0;}.有一个分数数列: 23581321,,,求出这个数列前20项之和 1235813#include
double sum(int n){ int i;double part = 0;
for(i = 1.0;i <= n;i++)
part +=(1.0 / i);return 2 * n-part;}
int main(void){ printf(“%.18fn”, sum(20));return 0;} 将从键盘输入的偶数写成两个素数之和。#include
int a,b,c,d;
scanf(“%d”,&a);
for(b=3;b<=a/2;b+=2)
{
for(c=2;c<=sqrt(b);c++)
if(b%c==0)break;
if(c>sqrt(b))
d=a-b;
else break;
for(c=2;c<=sqrt(d);c++)
if(d%c==0)break;
if(c>sqrt(d))
printf(“%d=%d+%dn”,a,b,d);
} } 1:5位跳水高手参加10米高台跳水决赛,有好事者让5人据实力预测比赛结果.A选手说:B第二,我第三B选手说:我第二,E第四;C选手说:我第一,D第二;D选手说:C最后,我第三;E选手说:我第四,A第一.决赛成绩公布之后,每位选手的预测都只说对了一半,即一对一错.请编程解出比赛的实际名次.
1.#include
#include
} for(i=1;i<=4;i++)
b=a[0];a[0]=a[3];a[3]=b;b=a[1];a[1]=a[2];a[2]=b;for(k=1;k<=4;k++)printf(“%d”,a[k-1]);{ } a[i-1]=(a[i-1]+5)%10;a[4-j]=b%10;b=b/10;}
2、编写程序,对输入两个正整数m和n,求出它们的最大公约数和最小公倍数 #include “stdio.h” #include “math.h” void main(){ int n,m,maxgy,mingb,i,p;printf(“please input n and m:”);scanf(“%d%d”,&n,&m);if(n>m){
p=n;n=m;m=p;/*m和n交换*/ } for(i=n;i>=1;i--)
if(m%i==0&&n%i==0)
break;maxgy=i;printf(“nmaxgy=%d mingb=%dn”,maxgy,m*n/maxgy);
}
2、编写程序,对输入两个正整数m和n,求出它们的最大公约数和最小公倍数
#include “stdio.h” #include “math.h” void main(){ int n,m,maxgy,mingb,t,p;printf(“please input n and m:”);scanf(“%d%d”,&n,&m);if(n>m){
p=n;n=m;m=p;/*m和n交换*/ } p=m*n;while(m%n!=0){
t=m%n;
m=n;
n=t;} maxgy=n;printf(“nmaxgy=%d mingb=%dn”,maxgy,p/maxgy);} #include “stdio.h” #include “math.h” void main(){ int n,m,maxgy,p;int maxgy1(int m,int n);printf(“please input n and m:”);scanf(“%d%d”,&n,&m);if(n>m){
p=n;n=m;m=p;/*m和n交换*/ } p=m*n;
maxgy=maxgy1(m,n);printf(“nmaxgy=%d mingb=%dn”,maxgy,p/maxgy);} int maxgy1(int m,int n){ if(n==0)
return m;else return maxgy1(n,m%n);}
3输入n判断n是否为素数 #include “stdio.h” void main(){ int n,i,flag;flag=1;printf(“please input n:”);scanf(“%d”,&n);for(i=2;i if(n%i==0) { flag=0; break; } if(flag==1) printf(“n%d is ssn”,n);else printf(“n%d is not ssn”,n);} #include “stdio.h”、求100以内的所有素数,并按10个一行进行打印。#include “math.h” void main(){ int n,i,flag,sum;sum=0;for(n=2;n<=100;n++){ flag=1; for(i=2;i<=sqrt(n);i++) if(n%i==0) { flag=0; break; } if(flag==1) { sum++; printf(“%5d”,n); if(sum%10==0) printf(“n”); } } printf(“n”);} 4、用递归求1到100的和 #include return 1;else return n+lj(n-1);} 累加法求1到100的和 #include sum=sum+i;printf(“result=%dn”,sum);} #include printf(“result=%dn”,sum);} #include if(i>100) break; sum=sum+i; i++;} printf(“result=%dn”,sum);} 求20到40以及60到80中所有能被3整除数的和 #include if(i%3==0) sum=sum+i; if(i==40) i=i+19;} //i%3==0&&i>=20&&i<=40||i>=60&&i<=80 printf(“result=%dn”,sum);} A+aa+aaa+….+a…..a #include t=t*10+a; sum=sum+t;} printf(“result=%ldn”,sum);} 1、请从键盘上输入年、月、日,求该年月日是该年的第多少天? #include int year,month,day,sum=0,i;int days(int,int);printf(“please input year month and day:”);scanf(“%d%d%d”,&year,&month,&day);for(i=1;i int day;switch(month){ case 1: case 3: case 5: case 7: case 8: case 10: case 12:day=31;break;case 4: case 6: case 9: case 11:day=30;break;case 2:if(year%4==0&&year%100!=0||year%400==0) } return day; day=29;else day=28;} 2、求3到1000内所有尾数为3的素数之和。#include } int ss(int n){ int flag=1,i;for(i=2;i<=sqrt(n);i++)if(n%i==0){ flag=0;break;int i,sum=0;int ss(int);for(i=3;i<=1000;i++)if(ss(i)==1) if(i%10==3) sum=sum+i;printf(“result=%dn”,sum); } } return flag; 3、从键盘上输入一个整数,将它拆成质因子的乘积。例如 18=2*3*3 #include } if(n%i==0){ } else i++;printf(“%d*”,i);n=n/i;printf(“b n”);} 4、从键盘上输入一串字符,统计字母、数字、空格和其它字符的个数。#include char ch;int c,d,s,o;c=d=s=o=0;while((ch=getchar())!=10){ if(ch>='a'&&ch<='z'||ch>='A'&&ch<='Z') c++;else if(ch>='0'&&ch<='9') d++;else if(ch==' ') s++;else o++;} printf(“c=%dnd=%dns=%dno=%dn”,c,d,s,o);} 5、从键盘上输入10个数,求它们的最大值。#include } int n,i,max;scanf(“%d”,&n);max=n;for(i=1;i<10;i++){ scanf(“%d”,&n);if(max max=n;} printf(“max=%dn”,max);3. 一个数恰好等于它的因子之和,这个数就称为“完数”。例如,6的因子为1,2,3而6=1+2+3,因此6是完数。编程找出求1000以内的所有完全数。#include “stdio.h” #include “math.h” void main(){ int n,s,i,k;for(n=1;n<=1000;n++){ s=0; for(i=1;i if(n%i==0) s=s+i; if(n==s) { printf(“%5d its factors is ”,n); for(k=1;k if(n%k==0) printf(“%d,”,k); printf(“b n”); } } printf(“n”);} .打印出杨辉三角形(要求打印出10行如下图)#include “stdio.h” #include “math.h” void main(){ long i,j,x,y,z,k;long jc(int);for(i=0;i<10;i++){ for(j=0;j<=i;j++) { x=jc(i); y=jc(j); z=jc(i-j); printf(“%4d”,x/(y*z)); } printf(“n”);} } long jc(int n){ long x=1,k;if(n==0) return 1;else for(k=1;k<=n;k++) x=x*k;return x;} #include “stdio.h” #include “math.h” void main(){ long i,j,x,y,z,k;for(i=0;i<10;i++){ for(j=0;j<=i;j++) { x=y=z=1; for(k=1;k<=i;k++) x=x*k; for(k=1;k<=j;k++) y=y*k; for(k=1;k<=i-j;k++) z=z*k; printf(“%4d”,x/(y*z)); } printf(“n”);} } 7.用*打印图形 #include “stdio.h” void main(){ int i,j,k,n;printf(“please input n:”);scanf(“%d”,&n);for(i=1;i<=n;i++)/**/ { for(j=1;j<=40-i;j++)/**/ printf(“ ”); for(k=1;k<=2*i-1;k++)/**/ printf(“*”); printf(“n”);} for(i=n-1;i>=1;i--)/**/ { for(j=1;j<=40-i;j++)/**/ printf(“ ”); for(k=1;k<=2*i-1;k++)/**/ } } printf(“*”);printf(“n”);