风雨钝剑的小窝

Tag:个人编写的小程序，原始发布位置在matlab盲去卷积图像复原程序

07年根据老研学论坛上一个网友贴出的文献所写，只具有基本功能，其中paddedsize是dipum函数，如有转载，请注明出处。

——by realasking

代码

 matlab 
function image=bldconv(g)
%Blind-deconv example
%Created by realasking
%Welcome to visit my blog: hi.baidu.com/realasking
%finish time:5-22-2007
 
a1=0.1;
a2=0.01;
 
PQ=paddedsize(size(g));
G=fft2(g,PQ(1),PQ(2));
 
[y x]=size(G);
htemp=ones(3);
h0=freqz2(htemp,PQ(1),PQ(2));
R=Rcreat(y,x);
 
H=h0;
for k=1:10 %计算IMG和psf
iMG=(conj(H).*G)./((conj(H).*H)+a1*R);
H=conj(iMG).*G./((conj(iMG).*iMG)+a2*R);
end
 
IMG=mat2gray(real(ifft2(iMG)));
image=IMG(1:size(g,1),1:size(g,2));
imwrite(image,'testfile.tif')
imshow(image,[]);
hold on;
 
end
 
%===========================
%计算R矩阵的函数Rcreat
%===========================
 
function R=Rcreat(y,x)
%R矩阵生成子函数
%by Realasking
%为bdeconv.m编制
 
i=1:y
j=1:x
 
RI=zeros(y,x);
RJ=zeros(y,x);
R=zeros(y,x);
 
for k=1:y %向量化代码生成R的矩阵
RI(k,i)=-2*cos(2*pi*i./y);
end
 
for k=1:x
RJ(j,k)=-2*cos(2*pi*j'./x);
end
 
R=4+RI+RJ;
 
end

Tag:自己写的一个程序，原始发布在：matlab微分方程循环求解的一种方法，以前matlab功能还不完善的时候有用，现在已经没有什么用处了。

一个取巧的做法，以前自己实现的一个范例，如有转载，请注明出处。

——by realasking

代码

 matlab 
function A=test1
%微分方程求解，由realasking创建
%未经允许请勿转载
%如果经过允许转载时，请保留此信息
%2007年5月13日
q=[5;5;5;4;6;9];
t0=0;tf=20;
x0=[0 0 0 0]
m=size(q)
for i=1:size(q)  
fid=fopen('a.m','w')  
fprintf(fid,'function xdot=a(t,x)\n');  
fprintf(fid,'q=%d;\n',q(i,1));  
fprintf(fid,'xdot=zeros(4,1);\n');  
fprintf(fid,'xdot(1)=(1/45.4)*(1500*(x(2)-x(1))+22000*(x(4)-x(3))+192000*(q-x(3)));\n') ;
fprintf(fid,'xdot(2)=(1/317.5)*(1500*(x(1)-x(2))+(x(3)-x(4))*22000);\n') ;
fprintf(fid,'xdot(3)=x(1);\n') ;
fprintf(fid,'xdot(4)=x(2);\n') ;
fclose(fid) ;
[t,x]=ode23('a',t0,tf,x0)
A(:,:,i)=x;
delete a.m  
end 

Tag:自己写的一个程序，原始发布地址在：matlab处理大学物理实验数据程序之最小二乘法

以前自己根据公式写的，批改实验报告方便，也是以前在百度发布过的老文，如有转载，请注明出处。

——by realasking

代码

 matlab 
function [a,b,r,ua,ub]=line_fit(C,D)
%最小二乘法拟合y=ax+b直线，并计算a、b的A类不确定度
%程序编制：realasking
 
%个人博客：http://hi.baidu.com/realasking
%输入变量:C、D，两个列向量，分别是x，y
%输出变量:a,b,r,ua,ub,r是残差
tmp=(sum(C)^2)-size(C,1)*sum(C.^2);
b=(sum(C.*D)*sum(C)-sum(D)*sum(C.^2))/tmp;
a=(sum(C)*sum(D)-size(C,1)*sum(C.*D))/tmp;
r=D-(a*C+b);
sy=sqrt(sum(r.^2)/(size(C,1)-2));
sa=sy/sqrt(sum((C-mean(C)).^2));
sb=sy*sqrt(((mean(C)^2)/(sum((C-mean(C)).^2))+(1/size(C,1))));
t=[1 12.7;
    2 4.3;
        3 3.18;
      4 2.78;
    5 2.57;
        6 2.45;
      7 2.36;
    8 2.31;
        9 2.26;
      10 2.23;
    15 2.13;
        20 2.09;
      21 1.96]
if(size(C,1)<=10)
    tt=t(size(C,1),2);
    elseif(size(C,1)<=15)
        tt=t(10,2)+((t(11,2)-t(10,2))*((size(C,1)-10)/(15-10)));
  elseif(size(C,1)<=20)
      tt=t(11,2)+((t(12,2)-t(11,2))*((size(C,1)-15)/(20-15)));
else
    tt=t(13,2);
    end
    ua=tt*sa;
    ub=tt*sb;

参考文献

吴平，大学物理实验教程，机械工业出版社