对一实测信号的处理，希望大家参与讨论

甜心宝贝

songzy41 发表于 2009-11-9 08:15

                               
登录/注册后可看大图

LZ在做 STFT时可能直接用：
tfrstft(yy)
可以设置为：

双对数坐标轴和直接幅值去对数是有区别的吧？

shajue

回复 2 # songzy41 的帖子

请问一个问题，为什么你测量的时候，采600秒，这么长的数？是为了做FFT的时候提高频率分辨率吗？
对于这种桥梁信号，你这一次测量的结果，与下一次测量的结果，应该有很大的差别吧？（特别是幅值方面，因为桥梁上的车不一样）。如果差别很大的话，你这种测量的意义是否不大？你又是怎样解决的呢？
不知道我说得对不对？请指导。

bcyangbc

shajue 发表于 2011-6-8 21:24

                               
登录/注册后可看大图

回复 2 # songzy41 的帖子

请问一个问题，为什么你测量的时候，采600秒，这么长的数？是为了做FFT的时候提 ...

不好意思，因现在未做这个方面了，很久没有登陆论坛！

你说的很对，取这么长的数据就是为了提高频率分辨率！若再去测量一次，由于荷载或外界条件不同，可能幅值会有所不同，但对于成桥，频率不会改变。

而且，对于桥梁结构来说，我们更关心频率，不是吗？
不知道我说得对不对？

bcyangbc

回复 14 # lianlian123 的帖子

用采样率为5Hz的GPS接收机测的，不存在问题，这方面有很多研究。

bcyangbc

回复 10 # MinnieLMickey 的帖子

程序已共享，不知是否看到？

bcyangbc

wildcranehust 发表于 2011-3-13 17:55

                               
登录/注册后可看大图

回复 35 # bcyangbc 的帖子

lz你好，哪个rawy.dat文件该用什么软件打开才能看见里面的数据结构那，我下载完 ...

rawy.dat是Matlab的数据文件，在matlab中直接load即可？
您没看代码？

lg12170226

请问楼主  我也是对物体振动测量 我滤波后做fft  频谱图有不止一个峰值  输入是单正弦信号，是不是我滤波有问题  询问一下Vondrak滤波法  谢谢

comoma2011

谢谢lz，对于一般的实测信号（必然包括噪声），都需要滤波处理么？

redplum

由于没有原程序啊

rwwttt

回复 1 # bcyangbc 的帖子

我觉得，关于实信号分析，相对来说还是中值滤波效果好一点，

sanmuhaifeng

回复 1 # bcyangbc 的帖子

1. N=length(y);
   
2. fidelity=1/epsilon;
   
3. P=ones(1,N);
   
4. a=zeros(1,N-3);
   
5. b=zeros(1,N-3);
   
6. c=zeros(1,N-3);
   
7. d=zeros(1,N-3);
   
8. for i=1:N-3
   
9.     a(i)=6*sqrt(t(i+2)-t(i+1))/((t(i)-t(i+1))*(t(i)-t(i+2))*(t(i)-t(i+3)));
   
10.     b(i)=6*sqrt(t(i+2)-t(i+1))/((t(i+1)-t(i))*(t(i+1)-t(i+2))*(t(i+1)-t(i+3)));
    
11.     c(i)=6*sqrt(t(i+2)-t(i+1))/((t(i+2)-t(i+1))*(t(i+2)-t(i))*(t(i+2)-t(i+3)));
    
12.     d(i)=6*sqrt(t(i+2)-t(i+1))/((t(i+3)-t(i+1))*(t(i+3)-t(i+2))*(t(i+3)-t(i)));
    
13. end
    
14. S=0;
    
15. for i=1:N-3
    
16.     S=S+(a(i)*y(i)+b(i)*y(i+1)+c(i)*y(i+2)+d(i)*y(i+3))^2;
    
17. end
    
18. % a
    
19. % b
    
20. % c
    
21. % d
    
22. AM=zeros(N,N);
    
23. for i=1:1:N
    
24.     for j=-3:1:3
    
25.         if (j+i<=0 || j+i>=N+1)
    
26.             %AM(i,j+i)=0;
    
27.             continue;
    
28.         end
    
29.         switch j
    
30.             case -3
    
31.                 if i-3<=0 || i-3>=N-2
    
32.                     ai3=0;
    
33.                     di3=0;
    
34.                 else
    
35.                     ai3=a(i-3);
    
36.                     di3=d(i-3);
    
37.                 end
    
38.                 AM(i,i+j)=ai3*di3;
    
39.             case -2
    
40.                 if i-2<=0 || i-2>=N-2
    
41.                     ai2=0;
    
42.                     ci2=0;
    
43.                 else
    
44.                     ai2=a(i-2);
    
45.                     ci2=c(i-2);
    
46.                 end
    
47.                 if i-3<=0 || i-3>=N-2
    
48.                     bi3=0;
    
49.                     di3=0;
    
50.                 else
    
51.                     bi3=b(i-3);
    
52.                     di3=d(i-3);
    
53.                 end
    
54.                 AM(i,i+j)=ai2*ci2+bi3*di3;
    
55.             case -1
    
56.                 if i-1<=0 || i-1>=N-2
    
57.                     ai1=0;
    
58.                     bi1=0;
    
59.                 else
    
60.                     ai1=a(i-1);
    
61.                     bi1=b(i-1);
    
62.                 end
    
63.                 if i-2<=0 || i-2>=N-2
    
64.                     bi2=0;
    
65.                     ci2=0;
    
66.                 else
    
67.                     bi2=b(i-2);
    
68.                     ci2=c(i-2);
    
69.                 end
    
70.                 if i-3<=0 || i-3>=N-2
    
71.                     ci3=0;
    
72.                     di3=0;
    
73.                 else
    
74.                     ci3=c(i-3);
    
75.                     di3=d(i-3);
    
76.                 end
    
77.                 AM(i,i+j)=ai1*bi1+bi2*ci2+ci3*di3;
    
78.             case 0
    
79.                 if i-3<=0 || i-3>=N-2
    
80.                     di3=0;
    
81.                 else
    
82.                     di3=d(i-3);
    
83.                 end
    
84.                 if i-2<=0 || i-2>=N-2
    
85.                     ci2=0;
    
86.                 else
    
87.                     ci2=c(i-2);
    
88.                 end
    
89.                 if i-1<=0 || i-1>=N-2
    
90.                     bi1=0;
    
91.                 else
    
92.                     bi1=b(i-1);
    
93.                 end
    
94.                 if i>=N-2
    
95.                     ai=0;
    
96.                 else
    
97.                     ai=a(i);
    
98.                 end
    
99.                 AM(i,i+j)=ai^2+bi1^2+ci2^2+di3^2+epsilon*P(i)/(N-3);
    
100.             case 1
     
101.                 if i-2<=0 || i-2>=N-2
     
102.                     ci2=0;
     
103.                     di2=0;
     
104.                 else
     
105.                     ci2=c(i-2);
     
106.                     di2=d(i-2);
     
107.                 end
     
108.                 if i-1<=0 || i-1>=N-2
     
109.                     bi1=0;
     
110.                     ci1=0;
     
111.                 else
     
112.                     bi1=b(i-1);
     
113.                     ci1=c(i-1);
     
114.                 end
     
115.                 if i>=N-2
     
116.                     ai=0;
     
117.                     bi=0;
     
118.                 else
     
119.                     ai=a(i);
     
120.                     bi=b(i);
     
121.                 end
     
122.                 AM(i,i+j)=ai*bi+bi1*ci1+ci2*di2;
     
123.             case 2
     
124.                 if i-1<=0 || i-1>=N-2
     
125.                     bi1=0;
     
126.                     di1=0;
     
127.                 else
     
128.                     bi1=b(i-1);
     
129.                     di1=d(i-1);
     
130.                 end
     
131.                 if i>=N-2
     
132.                     ai=0;
     
133.                     ci=0;
     
134.                 else
     
135.                     ai=a(i);
     
136.                     ci=c(i);
     
137.                 end
     
138.                 AM(i,i+j)=ai*ci+bi1*di1;
     
139.             case 3
     
140.                 if i>=N-2
     
141.                     ai=0;
     
142.                     di=0;
     
143.                 else
     
144.                     ai=a(i);
     
145.                     di=d(i);
     
146.                 end
     
147.                 AM(i,i+j)=ai*di;
     
148.         end
     
149.     end
     
150. end
     
151. AM=inv(sparse(AM));
     
152. By=zeros(1,N);
     
153. B=epsilon*P/(N-3);
     
154. for i=1:N
     
155.     By(i)=B(i)*y(i);
     
156. end
     
157. By=By';
     
158. fy=AM*By;
     
159. t_o=t;
     
160. y_o=fy;
     
161. figure(1);
     
162. subplot(211);
     
163. plot(t,yy,'b-');hold on;
     
164. plot(t_o,y_o,'r-','linewidth',2);
     
165. legend('原始数据','滤波结果','orientation','horizontal');
     
166. subplot(212);
     
167. plot(t,y_o-yy);
     
168. subplot(211);title('滤波');ylabel('位移(m)');
     
169. subplot(212);title('滤波残差');xlabel('时间(s)');ylabel('位移(m)');

复制代码

你好  请问一下，你在滤波去噪平滑的过程中，最长那段for语句 主要做了一个什么工作啊？  太长了，  不是很好理解   

sanmuhaifeng

回复 1 # bcyangbc 的帖子

谢谢大神，  本人是初学者，还有个问题想请教您，就是vondrak 滤波法在滤波时候在信号的截断频带上有较好的分辨率，  用着这边做信号的平滑处理好处体现在哪里咧？   会不会把高频的有用信号滤掉？另外有个问题我在考虑，在实际信号处理的过程中，到底是先去直流分量，趋势项后再滤波，还是先滤波再去直流分量？  

zhaowei_huang

支持一下 刚踏入这一行业

manguoyong

这个案例很好，学到了不少东西

manguoyong

很好，但是有点问题