从DTFT到DFT,计算频谱,并由频谱反求时间样点.docVIP

从DTFT到DFT，计算频谱，并由频谱反求时间样点 已知离散时间信号的离散时间傅里叶变换（DTFT）的表达式为： (DTFT正变换) 对频率离散化，在[0,2π]范围内取离散的频率样点，其中k=0,1,2,…,M-1，得 把记作X(k)，则有 (DFT正变换) 令，则上式可改写成： 则 经过矢量矩阵代换后，上式变为 X = x W (DFT正变换) 其中W是N行M列矩阵，且易证WWH=M IN。x是离散时间样点构成的1×N矢量，X是x的DTFT在一个周期上均匀取M个样点得到的离散频率样点矢量（DTFT是连续的）。 根据矩阵W中各行矢量间的正交性，易得反变换为: (DFT反变换) 于是，可以由离散的时间序列求得其DFT（可采用快速傅里叶变换FFT算法），再除以fs后取一个周期内的主值(k=[-M/2 M/2-1])便得到了x(t)的频谱；只不过自变量k对应的实际频率是k/(MTs) ＝k fs/M Hz。用matlab语言来描述，大致可写成下例的样子： Ts=0.0025; N=64; T=N*Ts; fs=1/Ts; %fs是抽样频率,必须满足抽样定理 M=N; %M是频率样点个数，这里简单取它跟时间样点数一样。 df = fs/M; %频谱分辨率 X=fft(x)/fs; %或X=x*W（须先计算出W来） f=(-M/2:M/2-1)*df; figure; plot(f, fftshift(abs(X))); %用fft得到的k=1,2,…,M对应离散DTFT值，用fftshift使右半边移到左边（利用DFT的周期性也是DTFT的周期性） title(幅频特性); 由于fft计算得到的是信号频谱乘以fs后的离散样值，要由ifft得到原始信号的离散样值，就必须对信号真实频谱乘以fs后再进行ifft运算。比如不能直接用： x = ifft(X); 这样来计算x。而必须使用： x=ifft(fs*X); %或x=fs*ifft(X); 这点很细节，看起来也很自然，但是必须牢记；且须注意到X是x的真实频谱。 之所以强调上面这点，是因为以下程序就用到了它： 其中Sussb是上边带已调信号，Hlow是对应相应基带消息信号的低通滤波器。①得到的确实是y2的频谱。②得到的DEM1也是真实的对相干解调信号进行频域低通滤波之后的真实信号频谱。根据上面讲到的，结果在③用ifft求对应时间信号时就必须先乘以fs。 这里看起来像是因为在①里乘了fs这里必须除以便进行补偿。其实不是这样。在之前设计低通滤波器Hlow时，假如给的是该滤波器的单位冲激响应h，用fft求频谱时就应该是：Hlow=fft(h)/fs;也除了一个fs,要是补偿的话，在③处岂不是应该乘以fs2来补偿？可以验证这样做是错误的。 y1=awgn(Sussb,20,measured); %调制信号通过AWGN信道,信噪比为20dB y2=y1.*cos(2*pi*fc*t); %相干解调 Y2=fft(y2)./fs; %①解调信号的频谱 DEM1=Y2.*Hlow; %②解调信号通过低通滤波器 dem1=real(ifft(DEM1*fs)); %③最终得到的解调信号 subplot(2,2,4) plot(t,dem1) title(信噪比为20dB时的解调信号)