Matlab之采样定理

haidaowang

$ d  P5 a- B9 F% n1 H当你看到我的博客的时候说明你已经做到快要吐血了，不过不用担心。可以看看我的，做的不一定对，不过可以借鉴一下！

    对连续信号进行等间隔采样形成采样信号，采样信号的频谱是原连续信号的频谱以采样频率为周期进行周期性的延拓形成的。设连续信号的的最高频率为Fmax，如果采样频率Fs>2Fmax，那么采样信号可以唯一的恢复出原连续信号，否则Fs<2Fmax会造成采样信号中的频谱混叠现象，不可能无失真地恢复原连续信号。
6 t7 c: Z# B/ d, p+ _& v/ r( @; Q
    首先是得到原信号的波形和频谱：
5 l- w$ o5 ?6 B
1 q: z7 a  H* t, ^6 M7 Lfunction Lab5_2_Sample_theorem1
! \$ T/ f! }7 a%t = 0:0.00001:0.1;
% {1 u/ H) X# ]%f = sin(2*pi*60*t)+cos(2*pi*25*t)+cos(2*pi*30*t);
  d! h8 [* f4 x2 N' R- @- X( E%plot(t,f,'r');
* x0 B% k8 ?# Q) b3 p: {fy = 'sin(2*pi*60*t) + cos(2*pi*25*t) +cos(2*pi*30*t)';
* K, |: t; C0 _# c  P. ?) E% 第一个输入变量函数是原信号函数，信号函数fy以字符串的格式输入
8 U9 Y6 K# f. a* s$ A9 x7 `6 d7 q% ~fs0 = 10000;
; Y' p3 Y: Y& Atp  = 0.1;
t   = -tp:1/fs0:tp;  001个样本点
3 i6 K- A. Q' D) C8 H5 X  s& {; n: Q+ jk1  = 0:999;
k2  = -999:-1;
m1  = length(k1);
; c6 F# k% i5 I( C4 c7 hm2  = length(k2);
, e" u( L  Q* _( pf   = [fs0*k2/m2,fs0*k1/m1];%f 1999
w   = [-2*pi*k2/m2,2*pi*k1/m1];%w = 2*pi*K/N,这里少算了头尾两点
) U; {  X3 P3 d+ P) s! Nfx1 = eval_r(fy);  %fx1 = sin(2*pi*60*t)+cos(2*pi*25*t)+cos(2*pi*30*t)
8 E" r: J* j9 ?7 H6 g) Z6 @9 \. DFX1 = fx1*exp(-1i*((1:length(fx1))./10)'*w); %除10好像得到了补0的效果
) D6 [: B) E) B%求原信号的离散时间傅里叶变换
figure(1);
, ?2 N' u5 _) n% I6 o6 Z5 J- X4 q%h画出原信号波形
! h+ ]9 H( K6 A1 z: w: gsubplot(1,2,1);plot(t,fx1,'r');
title('原信号');
- q" A+ `0 Y1 Y7 d: t; n8 Mxlabel('时间t/s');
axis([min(t),max(t),min(fx1),max(fx1)]);       %画原信号幅度频谱
subplot(1,2,2);
plot(f,abs(FX1),'r');
  e( u$ H5 Z& a& Atitle('原信号幅度频谱');
- Y2 ^8 r( x. J& ?; T% P: ?" txlabel('频率f');
5 V+ T; V7 B) w4 d# [& b0 l+ E7 V, Laxis([-1000,1000,0,max(abs(FX1)+5)]);          %画原信号幅度频谱

    其实把上面代码改一改就可以用到求采样信号的了。这里的卷积我试过很多种，以下的是我目前试出来唯一一个可以用的。

5 H7 [4 Y7 v5 e. a% J+ w% 对信号进行采样
fs   = 150;   % 采样频率
Ts   = 1/fs;
ws   = Ts*w;
ts   = Ts*(1:length(t));
fst  = sin(2*pi*60*ts)+cos(2*pi*25*ts)+cos(2*pi*30*ts);
Fs1 = fst*exp(-1i*((1:length(fst))./10)'*ws);
0 X, X2 |2 d6 S1 t5 Qlength(fst)
0 R3 a' w6 E' `2 [figure(2);
# [" O3 q7 Z( e  P, I' Qsubplot(1,2,1);plot(ts,fst,'r');
title('采样信号');
xlabel('时间t/s');
+ d% H  q6 R' o  Jaxis([min(ts),0.5,min(fst),max(fst)]);       %画原信号幅度频谱
subplot(1,2,2);
4 h/ H$ I$ J8 w/ x, D) C$ vplot(f,abs(Fs1),'r');
' T9 _, N. v3 l5 U5 h/ ptitle('幅度频谱');
' S- H9 p7 M% s! |  Txlabel('频率f');
. e. o. P* }1 wt1 = 1:length(t);%这里的冲击信号这么大是为了和采样信号的相匹配，最后结果因为用的是卷积所以这样
: X0 n! W8 |* k- G" C) |/ s5 s%求h(t)
ht = Ts*sin(2*pi*t1/Ts)./(t1*pi);
Ft = conv(ht,fst);
n = 1: max([length(fst)+length(ht)-1,length(fst),length(ht)]);
; x" u- s8 t' K5 d4 \figure(4);
plot(n,Ft,'b');
# R7 g$ G$ a/ maxis([200,450,min(Ft),max(Ft)]);
0 O& z) R* y) c# i# ~& ?( Q' s
    当然自己求卷积的话还有很多方法。一下就介绍一种方法：
: e: Z' T% A2 ~0 g: F" b% y" b3 d
y = -0.1:0.001:0.1;
n = -0.1/T:0.1/T;
* \' w4 b  N' V6 _$ ^M=ones(length(n),1)*y-n'*T*ones(1,length(y));
- R+ b% e( [7 L* ~/ z6 Efs = fg*sinc(frequency*M);
: Z4 ~7 s* C! L5 L* G/ @
    这个就是里面卷积常用的公式。就是有上下限，再相乘的，以上就是得到冲激函数。

yin123

Matlab之采样定理