werner 反褶积的程序

800522

我是第一次来论坛，请那位高手帮帮我，帮我找一些沃纳反褶积的资料和程序，谢谢！！！

风花雪月

个人只找到matlab和vc的，没发现fortran的

800522

那把你找到的matlab和vc的能否让我看看，谢谢阿

风花雪月

原帖由 800522 于 2008-1-15 11:23 发表

                               
登录/注册后可看大图

那把你找到的matlab和vc的能否让我看看，谢谢阿

留下你的email

800522

wangyani0704@yahoo.com.cn

风花雪月

1. function o= ls_inv_filter(w,wavelet,mu)
   
2. %LS_INV_FILTER: Spiking deconvolution
   
3. %
   
4. % This program computes an inverse filter (Spiking filter)
   
5. % from an estimate of a seismic wavelet. The actual
   
6. % output (convolution of the filter with the wavelet) is
   
7. % also computed.
   
8. %
   
9. % [f,o] = LS_INV_FILTER(w,NF,Lag,mu)
   
10. %
    
11. % IN w: the wavelet
    
12. % NF: lenght if the inverse filter
    
13. % Lag: the position of the spike in the desired output
    
14. % mu: Prewhitening in %
    
15. %
    
16. % OUT f: the filter
    
17. % o: the ouput or convolution of the filter with
    
18. % the wavelet
    
19. %
    
20. % Example:
    
21. %
    
22. % w = [2,1]; % the wavelet
    
23. % [f,o] = ls_inv_filter(w,20,1,2); % the filter and the output
    
24. % figure(1); plot(f);title('Filter')
    
25. % figure(2); plot(o);title('Wavelet \otimes Filter')
    
26. 
    
27. NW = max(size(w)); % lenght of the wavelet
    
28. 
    
29. NF=max(size(wavelet));
    
30. [mc,mr]=size(wavelet);
    
31. if mc <= mr; wavelet = wavelet'; end;
    
32. 
    
33. [mc,mr]=size(w);
    
34. if mc <= mr; w = w'; end;
    
35. 
    
36. b=conv(wavelet,w);
    
37. C = convmtx(w,NF); % Convolution matrix
    
38. 
    
39. R = C'*C+mu*eye(NF)/100.; % Toeplitz Matrix
    
40. rhs = C'*b; % Right hand side vector
    
41. f = inv(R)*rhs; % Filter
    
42. 
    
43. o = conv(f,w); % Actual output
    
44. % o=o(NF:NW+(NF-1),:);
    
45. % o=o(1:NW,:);
    
46. n=floor(NF/2);
    
47. o=o(n:NW+(n-1),:);  

复制代码