请问如何识别出这个信号里面的低功率部分的频率。

Laplacer

我也不知道我标题的表述对不对。。。。
问题是这样
x=1.5*sin(2*pi*151*t)+0.001*sin(2*pi*501*t)+r(t)
r(t)是一个正态随机信号。
我需要画出这个信号的功率谱密度。
我的代码如下：

1. 
   
2. clc; clear all;
   
3. N=2000;T=1;fs=N/T;df=fs/N;f=(0:N-1)*df;
   
4. t=linspace(0,T,N);
   
5. x=1.5*sin(2*pi*151*t)+0.001*sin(2*pi*501*t)+randn(1,length(t));
   
6. w1=hanning(N);sf2=N/sum(w1.^2);
   
7. x1=x.*w1';y1=fft(x1)/N;
   
8. p1=abs(y1).^2/df*2*sf2;
   
9. subplot(2,1,1);
   
10. semilogy(f(1:N/2),p1(1:N/2))
    
11. hn=hanning(N/8)';
    
12. noverlap=N/16;
    
13. [p2,f2]=pwelch(x,hn,noverlap,N/8,fs,'onesided');
    
14. subplot(2,1,2);semilogy(f2,p2)
    
15. 
    

复制代码

得到的曲线如下：

问题在于，x(t)里面有  0.001*sin(2*pi*501t)这个分量，我需要在图里面能看出来这个频率，就像能明显的看出 1.5*sin(2*pi*151t)一样。

求教啊。。。。

不是钱林

幅值太小。

Laplacer

不是钱林 发表于 2015-4-24 09:42
幅值太小。

嗯，就是幅值太小。。。我们老师说给我们出的难题。。。我也是醉了，这分量的功率加到信号上面，因为有噪声的影响，就是没区别。。。。我们boss非说能弄出来。。

不是钱林

滤波降噪，在那个频段上再细化谱分析一下。或是降噪后用小波分解或EMD分频段看一下，应该有

Laplacer

不是钱林 发表于 2015-4-24 10:02
滤波降噪，在那个频段上再细化谱分析一下。或是降噪后用小波分解或EMD分频段看一下，应该有

呃，我研究研究

Traly

幅值太小，不好出来啊，还有噪声

Laplacer

Traly 发表于 2015-4-24 15:08
幅值太小，不好出来啊，还有噪声

我们老师就是这么BT。。。

wanyeqing2003

我想楼主需要解决几个问题：
1、离散数据分析方法；
2、信号分析的误差；
3、数据量级与精度的关系。

我看楼主提供的信号可以分为三个部分：1、151频率成分；2、501频率成分；3、随机数据r。
其中的量级有三个：1、151频率的1.5；2、501频率的0.001；3、1.0的随机数据。

所以：
1、数据分析时，样本数据有限，会有截断误差。最大频谱10%以下的信号可能会被这个误差掩盖了，是不可靠的。
2、501频率成分要比151频率成分小1500倍，这样小的信号已被汹涌的误差掩盖无几，无法判别。
3、再加上随机信号的影响，501频率成分就更难找到。

个人见解，仅供参考。

wanyeqing2003

Laplacer 发表于 2015-4-24 18:21
我们老师就是这么BT。。。

这是一个理论上的数据，可否这样试一下。
先不考虑随机数，加大样本长度，例如样本总数取到10240或更大。

既然有要求，不妨做些尝试。大家都会的东西，再研究也就没多大意义了。呵呵。

Laplacer

wanyeqing2003 发表于 2015-4-26 09:59
这是一个理论上的数据，可否这样试一下。
先不考虑随机数，加大样本长度，例如样本总数取到10240或更大 ...

C:\Users\honeg\Desktop\N=10240.png
不加随机项，我算过了，因为y采用对数轴，所以看的出来。
然后N=10240的话，依然没有~~确实幅值太小，老师给的提示是“研究研究窗函数”我也是醉了。C:\Users\honeg\Desktop\N=2000无随机项.png

小逸

这。。。。。。。。。。。。

VibrationMaster

1.如果你们老师允许你使用任意长度的数据，那么可以用时域平均。也就是把数据按501那个成分的周期整数倍的长度，让后对截断后的数据平均，多次平均之后可以把噪声和151的成分去掉

hcharlie

VibrationMaster 发表于 2015-4-27 21:53
1.如果你们老师允许你使用任意长度的数据，那么可以用时域平均。也就是把数据按501那个成分的周期整数倍的 ...

取501整数周期不加窗（加矩形窗）是个好主意，多次平均可以抑制噪声，频域平均也一样。但如果要抑制到1/1000的水平，需要平均一百万次。希望你有耐心。你这是理想情况的数字模拟，真正工程问题平均1000000次也不一定能解决问题，因为时域平均要求频率绝对准确，工程上501的频率会有误差和不稳定，不容易取到完全等同相位的整数周期，另外噪声也不可能是理想正态分布，这还要经过工程实际的考验。

wanyeqing2003

Laplacer 发表于 2015-4-26 18:24
不加随机项，我算过了，因为y采用对数轴，所以看的出来。
然后N=10240的话，依然没有~~确实幅值太小， ...

看来研究有些进展了。
结合楼上几位的讨论，我们是否可以理出这样的思路。
1、不考虑随机项时，加大样本数量，如用到10240，在对数坐标下可以看出501Hz频率成分。
2、如果考虑随机项，可以用多段平均方法减少随机因素的影响。不过501Hz量级太小，可以试试。
3、采用加窗修正，目的是为了减小泄露的偏差。（我认为：加窗是为了美化曲线形态，是人为修正。加大样本数效果更好。）

VibrationMaster

1. 如果不想让151那个分量看到，再平均足够多次后，做谱图显示时把频率范围设成501附件即可
2.上面的做法你感觉掩耳盗铃。稍微高明一点的掩耳法是做低通滤波，把151的那个分量滤掉
3.为了在把501峰从噪声背景中突出来，不需要一百万次平均，因为FFT也相当于平均，所以如果是1024点FFT的话，大约用(1/0.001)^2/sqrt(1024)=3*10^4平均就差不多了
4.这种数据“操弄”不用加窗