自编程序分享-Wigner-Ville- Distribution

pengzk

1. function [WVD] = WignerVille(Sig,SampFreq,FreqBins);
   
2. %           Sig      : 输入信号
   
3. %    FreqBins       : 频率轴划分区间数(默认为512)
   
4. %  SampFreq       : 信号的采样频率
   
5. %         WVD       : 计算结果
   
6. 
   
7. if (nargin < 1),
   
8.     error('At least one input required');
   
9. end;
   
10. Sig = hilbert(real(Sig));        % 计算信号的解析信号
    
11. SigLen = length(Sig);           %  获取信号的长度   
    
12. if(nargin < 3),
    
13.     FreqBins = 512;
    
14. end
    
15. if(nargin < 2),
    
16.     SampFreq = 1;
    
17. end
    
18. FreqBins = 2^nextpow2(FreqBins);  
    
19. FreqBins = min(FreqBins, SigLen);
    
20. WVD = zeros (FreqBins,SigLen);                       
    
21. wait = waitbar(0,'Under calculation, please wait...');
    
22. for iIndex = 1:SigLen,
    
23.     waitbar(iIndex/SigLen,wait);
    
24.     MTau = min([iIndex-1,SigLen-iIndex,round(FreqBins/2)-1]);  
    
25.     tau = -MTau : MTau;
    
26.     Temp = rem(FreqBins + tau,FreqBins) + 1;
    
27.     WVD(Temp,iIndex) = Sig(iIndex + tau) .* conj(Sig(iIndex - tau));    %计算x*(t-1/2τ) x(t+1/2τ)
    
28. end;
    
29. close(wait);
    
30. WVD = fft(WVD)/FreqBins;
    
31. f = linspace(0,0.5,FreqBins)*SampFreq;
    
32. t = (0: SigLen-1)/SampFreq;
    
33. set(gcf,'Position',[20 100 400 300]);            
    
34. set(gcf,'Color','w');                                            
    
35. mesh(t,f,abs(WVD));
    
36. colormap('Jet');
    
37. colorbar;
    
38. axis([min(t) max(t) min(f) max(f)]);
    
39. ylabel('Frequency / Hz');
    
40. xlabel('Time / Sec');
    
41. title('Wigner-Ville distribution');

复制代码

pengzk

1. function TestWVD(),
   
2. 
   
3. % SampFreq = 500;
   
4. %
   
5. % t = 0:1/SampFreq:2;
   
6. %
   
7. % Sig = sin(50*2*pi*t+50*pi*t.^2);
   
8. SampFreq = 100;
   
9. t = 0 : 1/SampFreq : 15;
   
10. 
    
11. c1 = 2 * pi * 10;            % initial frequency of the chirp excitation
    
12. c2 = 2 * pi * 2.5;           % set the speed of frequency change be 1 Hz/second
    
13. c3 = 1 * pi * 1/1.5;
    
14. c4 = 2 * pi * -1/40;
    
15. 
    
16. Sig = sin(c1 * t + c2 * t.^2 / 2 + c3 * t.^3 /3 + c4 * t.^4 /4);   % get the A(t)
    
17. %Sig = Sig + sin((1.2*c1 * t + c2 * t.^2 / 2 + c3 * t.^3 /3 + + c4 * t.^4 /4));
    
18. 
    
19. 
    
20. figure(1)
    
21. plot(t,Sig);
    
22. axis tight
    
23. xlabel('Time/Sec')
    
24. 
    
25. figure(2)
    
26. Freq = linspace(0,SampFreq,length(t));
    
27. Spec = 2*fft(Sig)/length(t);
    
28. plot(Freq(1:end/2),abs(Spec(1:end/2)));
    
29. axis tight
    
30. xlabel('Frequency /Hz')
    
31. 
    
32. figure(3)
    
33. pause
    
34. WignerVille(Sig, SampFreq,512);

复制代码

geophysicszy

非常感谢！！！

lxb042

太帅了:handshake

asdf229955

非常感谢！非常感谢！:lol

nkdtxf

:victory: ，非常感谢啊，学习一下

robustcococole

学习了。。。。

shenlong1hao

非常谢谢！

zoezxc

Thank you for your sharing.

龙真诚

学习一下

hmlz

非常感谢楼主！！！！！

nianqingcuo

太厉害了！

gb9813

恩，学习一下

知识乞丐

非常感谢楼主！

bangbangni

很好，楼主很 牛