带限(窄带)信号傅里叶变换_依玲_百度空间

%--------------------------------------------------------------------------
% 带限(窄带)信号傅里叶变换
% ---------------
% 分析窄带信号的,带限谱数据变换到时间域问题.
% 关键词: FT,DFT,FFT,IFFT,窄带信号,傅立叶逆变换,扫频信号,频频率信号
% 目的: 把窄带的频域数据变换到时间域.使用DFT实现,或是IFFT实现的方法.
% 分析: 认识采样定理, 正负频关系,幅度关系.
% 意义: 认知离散FT变换一些细节问题.
%      理论的知识点虽然简单,但是在实现过程中会遇到许多问题.包括: dt,df,N,NF,N-1,NF-1,绘图等概念及技术.
% ---------------
% Coded by yiling
% Email:
% Department: College of Geo-Exploreration science and Technology of Jilin University, ChuangChun,China.
% Date: 2010-02-14
%--------------------------------------------------------------------------
close all;
clear;
clc;
%--------------
fs=5;   % GHz       采样频率
f0=0.5; % GHz       载频频率,即信号中心频率
t=0:1/fs:51;    % 原始信号时间轴
dt = 1/fs;      % 原始信号时间采样间隔
%
[tmp,N]=size(t);      % 原始信号时间采样点数
%-----------
FigWH = [283,137,480,320];    % 画图窗口大小
%--------------------------------------------------------------------------
% 载频信号
cs = cos(2*pi*f0*t);        % 余弦调制
%--------------------------------------------------------------------------
% 高斯信号 调制信号
tou = 20;
t0 = 0.75*tou;
u = exp(-4*pi*(t-t0).^2/tou^2);
u=u/max(u);
% 调制输出信号,即信号
s=u.*cs ;
%-----
% 绘制信号
figure('position',FigWH);
plot(t,u,'k--','LineWidth',2);
hold on;
plot(t,s,'r','LineWidth',2);
axis([t(1),t(end),-1.25,1.25]);
grid on;
ylabel('Normal Amplidute','FontSize',14); xlabel('Time (ns)','FontSize',14);
title('Origin Signal','FontSize',14);
set(gca,'FontSize',12);



%--------------------------------------------------------------------------
% 计算信号谱
NF=4096;                % 偶数, FFT变换长度
tr = [0:NF-1]*dt;       % 使用NF进行FT后,对应的时域来说,其对应的时间轴, dt与原信号相同,因为谱长度并没有改变(dt不变,谱长度即不变)
df=1/(dt*(NF-1));       % 频域 频率采样间隔,注意时间长度是dt*(NF-1),不是dt*NF. 但是差别不大,因为NF一般很大,NF与NF-1差别不大.
f=[-NF/2:NF/2-1]*df;    % 域域 频率轴,从负{zd0}频,到正{zd0}频, 可以证明,dt=1/(df*(NF-1)).
% 注意{zh0}采用这种方式形成频率轴,使用linspace(-fs/2,fs/2-df,NF),可能不会出现零频.
Uf=fftshift(fft(s,NF)); % 把零频放在中间
%----
figure('position',FigWH);
plot(f,abs(Uf)/max(abs(Uf)),'k','LineWidth',2);
axis([0,0.2*fs,-0.2,1.25]);
grid on;
ylabel('Normal Spectrum','FontSize',14); xlabel('Frequency (GHz)','FontSize',14);
title('Origin Spectrum','FontSize',14);
set(gca,'FontSize',12);



%--------------------------------------------------------------------------
%==========================================================================
%--------------------------------------------------------------------------
% 取带限信号
fid = find(0.3 < f & f < 0.7);      % 查找窄带谱有值的一段,正频段,在这里有值的一段是0.3GHz到0.7GHz
[tmp,Nsub]=size(fid);          
Nsub = ceil(Nsub/2)*2;              % 令取出的谱是2的整数倍,子带段的长度,正频
Nsub2= Nsub*2;                      % 子带段长度的两倍,正频
fid = fid(1):(fid(1)+Nsub-1);       % 取出的正频段的数编号
flmt = f(fid);                      % 取出谱的,频率轴
Ulmt = Uf(fid);                     % 取出的谱,正频
%-------------
figure('position',FigWH);
plot(flmt,abs(Ulmt)/max(abs(Ulmt)),'k','LineWidth',2);
axis([-inf,inf,-0.2,1.25]);
grid on;
ylabel('Normal Spectrum','FontSize',14); xlabel('Frequency (GHz)','FontSize',14);
title('SubBand positive Spectrum','FontSize',14);
set(gca,'FontSize',12);



%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% 重新补零但只是正频
Ulmtr = zeros(1,NF);
Ulmtr(fid) = Ulmt;                  % 谱长度出原谱一样,但只保留正频有谱值的一段
%---
figure('position',FigWH);
plot(f,abs(Ulmtr)/max(abs(Ulmtr)),'k','LineWidth',2);
axis([-inf,inf,-0.2,1.25]);
grid on;
ylabel('Normal Spectrum','FontSize',14); xlabel('Frequency (GHz)','FontSize',14);
title('ReBuild Total Spectrum(Only Positive)','FontSize',14);
set(gca,'FontSize',12);



%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% [1] 对重建谱进行逆变换
sr=ifft(ifftshift(Ulmtr));              % 谱长度出原谱一样,但只保留正频有谱值的一段,进行IFFT
% 由于只有正频谱,能量降低为一半,幅度降低为一半
%---
figure('position',FigWH);
plot(t,u,'k--','LineWidth',2);          % 绘制调制信号,即包络信号
hold on;
plot(t,u/2,'b--','LineWidth',2);        % 包络信号降低一半.
plot(tr,real(sr),'r','LineWidth',2);    % 绘制只用正频逆变换出的信号, % 取实部,
axis([tr(1),tr(N),-1.25,1.25]);
grid on;
ylabel('Normal Amplidute','FontSize',14); xlabel('Time (ns)','FontSize',14);
title('[1]ReBuild Signal from Positive Spectrum','FontSize',14);
set(gca,'FontSize',12);
legend('Modulate','Half Modulate','Reconstruct');



%--------------------------------------------------------------------------
% 对带限谱直接进行逆变换
fb = f(fid(1));     % 子带的起始频率,正频
fe = f(fid(end));   % 子带的终止频率,正频
fbe = f(fid);       % 子带频率轴
fbnormal = fb/(fs); % 子带的起始频率,正频,使用采样频率归一化
fenormal = fe/(fs); % 子带的终止频率,正频,使用采样频率归一化
fbenormal = linspace(fb,fe,Nsub);   % 子带频率轴,使用采样频率归一化
Ulmts = [conj(Ulmt(end:-1:1)),0,Ulmt]; % 子带谱,含正负频
%--
figure('position',FigWH);
plot([-fbe(end:-1:1),0,fbe],abs(Ulmts)/max(abs(Ulmts)),'k.','LineWidth',2);
axis([-inf,inf,-0.2,1.25]);
grid on;
ylabel('Normal Spectrum','FontSize',14); xlabel('Frequency (GHz)','FontSize',14);
title('ReBuild Sub Spectrum','FontSize',14);
set(gca,'FontSize',12);



%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% [2] DFT 实现, 使用DFT定义式进行反变换,频率轴:0->pi->2pi
% 时间长度为原始信号的长度N.即t,也可以为FT变换长度NF,即tr
% for ii=1:N
%     ss(ii) = sum(ifftshift(Uf).*exp(j*2*pi*[0:NF-1]/NF*(ii-1)));
% end
% ss = ss/NF; % 除以正变换的点数
%--------------------------------------------------------------------------
% [3] DFT 实现2, 使用DFT定义式进行反变换,频率轴:-pi->0->pi
% 时间长度为原始信号的长度N.也可以为FT变换长度NF
% for ii=1:N
%     ss(ii) = sum(Uf.*exp(j*2*pi*[-NF/2:NF/2-1]/NF*(ii-1)));
% end
% ss = ss/NF;
%--------------------------------------------------------------------------
% [4] DFT 实现, 使用DFT定义式进行反变换,代入频率轴进行计算,这与[3]计算方式相同.
% 时间长度为原始信号的长度N.也可以为FT变换长度NF
% for ii=1:N
%     ss(ii) = sum(Uf.*exp(j*2*pi*f/fs*(ii-1)));
% end
% ss = ss/NF;
%--------------------------------------------------------------------------
% [5] DFT 实现, 使用DFT定义式进行反变换,代入频率轴及时间轴进行计算,这与[3][4]计算方式相同.
% 时间长度为原始信号的长度N.也可以为FT变换长度NF
% for ii=1:N
%     ss(ii) = sum(Uf.*exp(j*2*pi*f/fs*(t(ii)/dt)));
% end
%--------------------------------------------------------------------------
% [6] DFT 实现, 使用DFT定义式进行反变换,代入频率轴及时间轴进行计算,这与[3][4]计算方式相同.
% 在时间域中采用新的时间采样间隔,形式新的时间轴tt,
tt = [0:2*N-1]*dt/2; % 新的时间轴
for ii=1:2*N
    ss(ii) = sum(Ulmt.*exp(j*2*pi*fbe*tt(ii)));                       % {dy}种实现,只对正频进行计算
%     ss(ii) = sum(conj(Ulmt).*exp(-j*2*pi*fbe*tt(ii)));                % {dy}种实现,只对负频进行计算. 与{dy}种实现等价(相等)
%     ss(ii) = ss(ii)*2;                                                % 由于只对正或负频进行变换,则计算结果幅度应加倍.
    ss(ii) = ss(ii) + sum(conj(Ulmt).*exp(-j*2*pi*fbe*tt(ii)));        % 叠加 负频的变换
end
ss = ss/NF;
% 结论: 对于实物理可实现信号,对其谱值的逆变换,只对正频变换,与只对负频变换是等价(相等的,变换后都是取其实部.real(ss))
% {zh1}幅度应该除以正变换的点数NF.
%--------------------------------------
% 绘图
figure('position',FigWH);
plot(t,u,'k--','LineWidth',2);
hold on;
plot(t,s,'b.','LineWidth',1);
plot(tt,real(ss),'r','LineWidth',2);        % 时间轴,如果是[2-5],时间轴为t(N点),如果是[1]时间轴为tr(NF点),dt相同. 如果是[6]时间轴为tt(),新dt同. % 取实部,
axis([t(1),50,-1.2,1.2]);
grid on;
ylabel('Normal Amplidute','FontSize',14); xlabel('Time (ns)','FontSize',14);
title('[6]ReBuild Signal from SubBand by DFT','FontSize',14);
set(gca,'FontSize',12);
legend('Modulate','Origin','Reconstruct');



%--------------------------------------------------------------------------
%==========================================================================
%--------------------------------------------------------------------------
% [7] 使用FFT,谱搬移到低频.谱长度不变(正负频共NF点,频率采样间隔为df).
% 则时间轴为tr(NF点),因为谱长度不变,则dt不变.
Uf = Uf*0;
Uf(NF/2+1:NF/2+Nsub) = Ulmt;        % 只有正频搬移到低频
%------
% 绘制搬移的谱
figure('position',FigWH);
plot(f,abs(Uf)/max(abs(Uf)),'k','LineWidth',2);
axis([0,0.2*fs,-0.2,1.25]);
grid on;
ylabel('Normal Spectrum','FontSize',14); xlabel('Frequency (GHz)','FontSize',14);
title('Spectrum After Migration to Low Position','FontSize',14);
set(gca,'FontSize',12);



%-------------------
% [7_1] 使用IFFT逆变换 对搬移到低频谱进行变换
ss = ifft(ifftshift(Uf));                % 对频谱搬移到低频段后(谱长度不变,则dt不变),直接采用IFFT变换,后再乘以移位项(下一行)
% 采用移位法,理论上应对正负频分析单独IFFT,并使用不同的移项,即正负频移位方向是不同.{zh1}再把两个移位后的时间信号相加.
% 这里负频谱已置零,则不需再IFFT及移位.
ss = ss.*exp(j*2*pi*fb*tr);            % 移位项.    
% 幅度不归一化,因为FFT与IFFT是对应的,两个变换长度相同, 由于只有正频,则幅度加倍,以与原信号对比.
ss = ss;
ss = ss*2;          % 由于只有正频谱,幅度需加倍
%-------
% 绘制变换结果
figure('position',FigWH);
plot(t,u,'k--','LineWidth',2);
hold on;
plot(t,s,'b.','LineWidth',1);
plot(tr,real(ss),'r','LineWidth',2);        % 取实部,
axis([t(1),50,-1.2,1.2]);
grid on;
ylabel('Normal Amplidute','FontSize',14); xlabel('Time (ns)','FontSize',14);
title('[7\_1]ReBuild Signal from LowShift SubBand by IFFT','FontSize',14);
set(gca,'FontSize',12);
legend('Modulate','Origin','Reconstruct');
text(1,-1.1,'Using IFFT by Full Lenght Spectrum','FontSize',14);



%--------------------------------------------------------------------------
% [7_2] 采用IFFT,谱搬移到低频,(df不变,共Nsub2+1个点),谱长度变短(正负频共Nsub2点,再加上一个零频点,共Nsub2+1个频点,频率采样间隔为df).
ss = ifft(ifftshift([conj(Ulmt(end:-1:1))*0,0,Ulmt]));
% 由于谱长度变短,则时间分辨率变低(dt变大).
dtnew = 1/(Nsub2*df);       % 或 dtnew = (NF-1)*df/(Nsub2*df)/((NF-1)*df) = dt*(NF-1)/Nsub2;
ts = [0:Nsub2]*dtnew;       % 新的时间轴,Nsub2+1个点
ss = ss.*exp(j*2*pi*fb*ts);      % 移位项,
% 幅度归一化为: (Nsub2+1)/NF, 由于只有正频,则幅度加倍,以与原信号对比.
ss = ss/NF*(Nsub2+1)*2;
% ****
% 上述[7_2]方法问题,由于谱长度变短,时间分辨率为小(dt变大),则ss无法表达高频信息,移位项在大dt下也无法表达高频信息,
% 即搬移到高频,但是采样率不足,表现不出高频信息.
%----------------------
figure('position',FigWH);
plot(t,u,'k--','LineWidth',2);
hold on;
plot(t,s,'b.','LineWidth',1);
plot(ts,real(ss),'r','LineWidth',2);
axis([t(1),50,-1.2,1.2]);
grid on;
ylabel('Normal Amplidute','FontSize',14); xlabel('Time (ns)','FontSize',14);
title('[7\_2]ReBuild Signal from SubBand by IFFT','FontSize',14);
set(gca,'FontSize',12);
legend('Modulate','Origin','Reconstruct');
text(1,-1.1,'Can not describe High Frequencn under big dt','FontSize',14);



%--------------------------------------------------------------------------
% [7_3] 使用IFFT,由于谱长度变短,(df不变,共Nsub2+1个点),时间分辨率为小(dt变大),则无法表达高频信息
% 则在IFFT后,先进行插值
ss = ifft(ifftshift([conj(Ulmt(end:-1:1))*0,0,Ulmt]));
% 由于谱长度变短,则时间分辨率变低(dt变大).
dtnew = 1/(Nsub2*df);       % 或 dtnew = (NF-1)*df/(Nsub2*df)/((NF-1)*df) = dt*(NF-1)/Nsub2;
ts = [0:Nsub2]*dtnew;       % 新的时间轴,Nsub2+1个点
% 用高时间分辨率时间轴tr(即FFT变换时间轴,NF个点)进行插值
sy = interp1(ts,ss,tr,'spline');        % 样条插值
sy = sy.*exp(j*2*pi*fb*tr);             % 移位项                  
% 因为FFT与IFFT长度不同,幅度归一化为: (Nsub2+1)/NF, 由于只有正频,则幅度加倍,以与原信号对比. 幅度归一化只与FFT时的变换点数及IFFT是变换点数有关.
sy = sy/NF*(Nsub2+1)*2;
%-----------
figure('position',FigWH);
plot(t,u,'k--','LineWidth',2);
hold on;
plot(t,s,'b.','LineWidth',1);
plot(tr,real(sy),'r','LineWidth',2);
axis([t(1),50,-1.2,1.2]);
grid on;
ylabel('Normal Amplidute','FontSize',14); xlabel('Time (ns)','FontSize',14);
title('[7\_3]ReBuild Signal from SubBand by IFFT','FontSize',14);
set(gca,'FontSize',12);
legend('Modulate','Origin','Reconstruct');
text(1,-1.1,'Using IFFT with bit-shift method','FontSize',14);
% 精度明显不比直接用DFT好,是否是插值问题??



%--------------------------------------------------------------------------
%--------------------------------- EOF -----------------------------------
%--------------------------------------------------------------------------



郑重声明：资讯 【带限(窄带)信号傅里叶变换_依玲_百度空间】由 发布，版权归原作者及其所在单位，其原创性以及文中陈述文字和内容未经(企业库qiyeku.com)证实，请读者仅作参考，并请自行核实相关内容。若本文有侵犯到您的版权， 请你提供相关证明及申请并与我们联系（qiyeku # qq.com）或【在线投诉】，我们审核后将会尽快处理。
—— 相关资讯 ——
鲁贯通电蒸汽硫化罐 电加热水硫化罐
鲁贯通电蒸汽硫化罐 电加热水硫化罐
2024/8/19
自助碾米机的优势用途
自助碾米机的优势用途
2019/11/7
广州腾洋造浪设备的基本原理及步骤
广州腾洋造浪设备的基本原理及步骤
2019/9/28
SSI编码器_{jd1}值编码器编程厂家_桁萱自动化科技(图)
SSI编码器_{jd1}值编码器编程厂家_桁萱自动化科技(图)
2019/5/11
郑州热交换器清洗除垢公司|梅州热交换器清洗除垢|上海神通化工清洗
郑州热交换器清洗除垢公司|梅州热交换器清洗除垢|上海神通化工清洗
2018/8/18
上海回收海力士没打字内存容量64GB
上海回收海力士没打字内存容量64GB
2023/3/25
3年期央票隔近两年重启,进一步增强紧缩预期- 浦东铜俱乐部—铜交易者的 ...
3年期央票隔近两年重启,进一步增强紧缩预期- 浦东铜俱乐部—铜交易者的 ...
2010/4/12
[转载]肉炒平菇做法+糖醋辣白菜，菊香茄子，大麦米甜豆粥<br>_西岭雪_ ...
[转载]肉炒平菇做法+糖醋辣白菜，菊香茄子，大麦米甜豆粥<br>_西岭雪_ ...
2010/4/21
冬天做激光祛斑好吗
冬天做激光祛斑好吗
2015/3/8
欧韵艺墅商铺售楼处电话【40推荐】
欧韵艺墅商铺售楼处电话【40推荐】
2015/3/16
总是脱发怎么办
总是脱发怎么办
2015/4/21
伊利牛奶价格如何，精品伊利牛奶批发【河南】
伊利牛奶价格如何，精品伊利牛奶批发【河南】
2016/5/25
SKF轴承的研发之路_唯⒈のsē綵.的空间_百度空间
SKF轴承的研发之路_唯⒈のsē綵.的空间_百度空间
2010/6/5
【掩面】= =+谢他娘的超生……_Pu言Pu语_百度空间
【掩面】= =+谢他娘的超生……_Pu言Pu语_百度空间
2010/5/16