作者：FOU
前言
本文在 http://blog.csdn.net/dznlong/article/details/2261150 的基础上，重新整理，并添加了具体的Matlab例子。最后进一步扩充到调制与解调。
（一）时域-频域变换的本质
我们知道，一个单波，由频率w，振幅A，相位b组成，在时域上（横轴表示时间，数轴表示振幅），可以用余弦（正弦）表示为：y=A*cos(2*pi*w*t+b)。
来看一个简单的例子。频率为50Hz的余弦波：
Fs = 1000; %采样频率，每秒采样1000点
T = 1/Fs; %每一点就代表1/1000秒
L = 1000; %一共采样1000个点（就是采样一整秒啦）
t = (0:L-1)*T; %每一点对应的时刻（从最开始的0秒，到最后的1秒）
x1 = cos(2*pi*50*t) ; %50Hz信号及采样
figure;plot(t,x1);title('x1 50Hz');
反过来，已知时域上的信号样子，怎样知道它是由哪（几）个频率的波组成的呢？
凭什么说一个任意信号一定能拆分几个波呢？因为傅里叶在1807年给出了证明，一个信号，无论多么复杂，都可以表示成几个余弦波的叠加（求和）。具体证明请维基百科。反正就是能拆分成各种（频率，振幅，相位）的波。
对上面这个图，一眼看上去，它就是一个单波，振幅是1，频率是。。。至少可以数一下，50个峰所以是50Hz。可是如果是50000Hz，难道我们要数到天亮吗？或者图形复杂一点，说不上来究竟是不是波峰，又如何计数呢？
从数学或者程序的角度，我们有一个很笨但是管用的方法，就是我构造各种频率，各种相位的余弦波，与信号点积求和，波形一致的得到的和最大（正值最大峰），不相关的接近0（0），相反的最小（负值最大峰）。
(a)确定频率w
先不考虑相位。
1. 把一个30Hz的余弦波，和原始50Hz余弦波各点相乘，结果在[-1 1]区间，直观的感觉，每出现一个1，就会有一个对应的-1和它抵消掉，最终求和为（接近）0。
aCos = cos(2*pi*30*t);
vCos = aCos.*x1;
sum(vCos)
figure;subplot(3,1,1);plot(x1);subplot(3,1,2);plot(aCos);subplot(3,1,3);plot(vCos);title('check 30Hz');
>> 1.8519e-13
aCos = cos(2*pi*50*t);
vCos = aCos.*x1;
sum(vCos)
figure;subplot(3,1,1);plot(x1);subplot(3,1,2);plot(aCos);subplot(3,1,3);plot(vCos);title('check 50Hz');
>> 500
3. 极端情况，如果匹配波频率是0呢？那它就是一条直线，振幅为1。结果就是原始波的各点简单相加。
aCos = cos(2*pi*0*t);
vCos = aCos.*x1;
sum(vCos)
figure;subplot(3,1,1);plot(x1);subplot(3,1,2);plot(aCos);subplot(3,1,3);plot(vCos);title('check 0Hz');
>> 1.7808e-13
来做一个全范围测试，写个循环。我们来构建从0Hz到500Hz的501个波（以1Hz为单位逐渐增加）。为什么只到500Hz呢？因为原始信号采样频率是1000Hz（每秒采1000个点），所以我们新构建的波，为了点积求和，必须保持长度一致，一秒也是1000个点；一个周期至少要有一正一负，1000个点最多只能表示到500Hz，再往上就不具有周期性，没意义了。
cs = []; %贮存结果
for i =0:500 %创建不同501个频率的cos波
aCos = cos(2*pi*i*t); %创建当前频率的波
vCos = aCos.*x1; %点积
cs = [cs sum(vCos)]; %求和
end
figure;plot(0:1:500,cs);title('cos 0Hz-500Hz');
我们看到，只有50Hz的余弦波，与原始的50Hz余弦波点积求和，有显著的大值，表现为正相关；其它Hz的余弦波，点积求和接近于0，不能得到有意义的结果，表现为不相关。为什么只有同样频率的波点积求和有响应呢？
cos(2*pi*w1*t).*cos(2*pi*w2*t) = (cos(2*pi*(w1+w2)*t)+cos(2*pi*(w1-w2)*t))/2;第一部分始终是周期波，具有对称性，正负各点互相抵消，点积求和为0。第二部分呢，绝大部分情况下，也是一个周期波，点积求和为0。当且仅当w1==w2时，频率为0，值为1，才可以积累出一个大值。
（b）确定振幅A
正相关显示求和结果为500，为什么是500呢？这个500又怎样对应到原始振幅上呢？
50Hz与50Hz点积相乘，得到的结果为50Hz波的平方：y = cos(2*pi*50*t).*cos(2*pi*50*t)=cos(2*pi*50*t).^2; 其中t=0:0.001:1。
sum(y) = sum(cos(2*pi*50*0/1000)^2+cos(2*pi*50*1/1000)^2+cos(2*pi*50*2/1000)^2+...+cos(2*pi*50*1000/1000)^2;
又因为sin(a)^2+cos(a)^2 = 1 => cos(pi/2-a)^2+cos(a)^2=1;所以对上面的式子中的每一项，总能找到另一项和它相加等于1。一共有1000个项，正好凑出500个对子，每个对子为1，所有求和为500。只要是在整周期上均匀分布的点，总能互相找到配对求和为0。不是均匀分布就悬了。
直接看上面推导公式，每个采样点供应了1/2，一共N=1000个采样点，所以最终求和为500。
更一般化，求和得到的值，在最理想（原始波和测试波频率、振幅都一致）的状况下，得到的求和结果就为 A*N/2，A为原始波的振幅，N为采样点。并且推出，采样点频率越高，或者采样时间越长，都会增加样点数从而增高最终的和。
从最终的和反推出原始振幅就很简单了，s=A*N/2 => A=s*2/N;
我们可以试一下，求和值只和样点数与原始振幅相关，和具体的频率无关。
比如30Hz点乘30Hz，点积求和为500：
x1 = cos(2*pi*30*t) ;
aCos = cos(2*pi*30*t);
vCos = aCos.*x1;
sum(vCos)
>> 500
比如300Hz点乘300Hz，点积求和还是500：
x1 = cos(2*pi*300*t) ;
aCos = cos(2*pi*300*t);
vCos = aCos.*x1;
sum(vCos)
>> 500
但是，但是，在两个极端情况，频率为极小值0Hz，或者频率为极大值500Hz时，点积求和却为1000，这又是为什么呢？
x1 = cos(2*pi*0*t) ;
aCos = cos(2*pi*0*t);
vCos = aCos.*x1;
sum(vCos)
>> 1000
x1 = cos(2*pi*500*t) ;
aCos = cos(2*pi*500*t);
vCos = aCos.*x1;
sum(vCos)
>> 1000
这是因为，当极小或极大时，点积之后的波形都变成了水平波，每个样点自己就都是1，不用再凑对子了。多少个样点就是多少和。0Hz的情况下，cos(0)^2 = 1^2 = 1，始终等于1；500Hz的情况下，因为这是当前采样频率下（1000个采样点/秒）能表示的最高频率波，波形就是[1,-1,1,-1....,1,-1]，它的平方也始终是1。在这两个极端情况下，s=A*N => A = s/N;
还有一个小问题，我们的测试波是从0Hz到500Hz，每次增加1Hz。聪明的读者会问，为什么每次增加1Hz呢？可不可以只增加0.1Hz，0.01Hz，这样会测的更准确些吗？我们来试试看。
cs = [];
for i =0:0.1:500 %步长为0.1
aCos = cos(2*pi*i*t);
vCos = aCos.*x1;
cs = [cs sum(vCos)];
end
figure;plot(0:0.1:500,cs);title('cos 0Hz-500Hz');
我们发现，在整数频率上不相关，但是再细分到小数频率上，表现出了各种不可忽略的相关性，这是怎么回事呢？
这是因为，如果测试波的频率不为整数的话，最后总会有个尾巴不是完整周期，其对称性被破坏，无法求和为0。所以测试波必须要有完整的周期。
那么怎么样能把测试精度提升到小数位置呢？我们的例子里，测试的是整一秒的采样，如果我们继续多采样几秒，比如采样2秒钟，那就是2000个采样点，那么所有0.5Hz（例如0.5,1.5,2.5)在一个周期有尾巴，在两个周期正好凑成一个整尾巴，又是一个整周期啦。所以2秒采样时，我们就精确到0.5Hz；3秒采样时，精确到1/3Hz。依此类推，精度=1/采样秒数。
上面的例子，我们想精确到0.1Hz，就要采样10整秒。
Fs = 1000; %采样频率，每秒采样1000点
T = 1/Fs; %每一点就代表1/1000秒
L = 10000; %一共采样10000个点（就是采样10整秒啦）
t = (0:L-1)*T; %每一点对应的时刻（从最开始的0秒，到最后的10秒）
x1 = cos(2*pi*50*t) ; %50Hz信号及采样
figure;plot(t,x1);title('x1 50Hz');
cs = [];
for i =0:0.1:500
aCos = cos(2*pi*i*t);
vCos = aCos.*x1;
cs = [cs sum(vCos)];
end
figure;plot(0:0.1:500,cs);title('cos 0Hz-500Hz 10 seconds');
看，所有的0.1Hz确实都正确表现出了不相关性。而且说明，频率的测试可以精确到任何精度，50Hz和50.1Hz，在1Hz精度下可能相差不大，但是采样时间一长，到0.1Hz下，它们就是彻底两个不同的波了。（注意，因为采样点多了10倍，50Hz正相关的求和值也相应增大了10倍。）
（c）确定相位b
最后，我们来给原始信号加上相位，然后看看怎样测出这个相位。
Fs = 1000; %采样频率，每秒采样1000点
T = 1/Fs; %每一点就代表1/1000秒
L = 1000; %一共采样1000个点（就是采样一整秒啦）
t = (0:L-1)*T; %每一点对应的时刻（从最开始的0秒，到最后的1秒）
x1 = cos(2*pi*50*t+pi/3) ; %50Hz信号及采样,由于相位关系，整体向左位移pi/3
figure;plot(t,x1);title('x1 50Hz pi/3');
前言
本文在 http://blog.csdn.net/dznlong/article/details/2261150 的基础上，重新整理，并添加了具体的Matlab例子。最后进一步扩充到调制与解调。
（一）时域-频域变换的本质
我们知道，一个单波，由频率w，振幅A，相位b组成，在时域上（横轴表示时间，数轴表示振幅），可以用余弦（正弦）表示为：y=A*cos(2*pi*w*t+b)。
来看一个简单的例子。频率为50Hz的余弦波：
Fs = 1000; %采样频率，每秒采样1000点
T = 1/Fs; %每一点就代表1/1000秒
L = 1000; %一共采样1000个点（就是采样一整秒啦）
t = (0:L-1)*T; %每一点对应的时刻（从最开始的0秒，到最后的1秒）
x1 = cos(2*pi*50*t) ; %50Hz信号及采样
figure;plot(t,x1);title('x1 50Hz');
反过来，已知时域上的信号样子，怎样知道它是由哪（几）个频率的波组成的呢？
凭什么说一个任意信号一定能拆分几个波呢？因为傅里叶在1807年给出了证明，一个信号，无论多么复杂，都可以表示成几个余弦波的叠加（求和）。具体证明请维基百科。反正就是能拆分成各种（频率，振幅，相位）的波。
对上面这个图，一眼看上去，它就是一个单波，振幅是1，频率是。。。至少可以数一下，50个峰所以是50Hz。可是如果是50000Hz，难道我们要数到天亮吗？或者图形复杂一点，说不上来究竟是不是波峰，又如何计数呢？
从数学或者程序的角度，我们有一个很笨但是管用的方法，就是我构造各种频率，各种相位的余弦波，与信号点积求和，波形一致的得到的和最大（正值最大峰），不相关的接近0（0），相反的最小（负值最大峰）。
(a)确定频率w
先不考虑相位。
1. 把一个30Hz的余弦波，和原始50Hz余弦波各点相乘，结果在[-1 1]区间，直观的感觉，每出现一个1，就会有一个对应的-1和它抵消掉，最终求和为（接近）0。
aCos = cos(2*pi*30*t);
vCos = aCos.*x1;
sum(vCos)
figure;subplot(3,1,1);plot(x1);subplot(3,1,2);plot(aCos);subplot(3,1,3);plot(vCos);title('check 30Hz');
>> 1.8519e-13
2. 把一个50Hz的余弦波，和原始50Hz余弦波各点相乘，结果在[0 1]区间。每个点的积都大于0，最终求和为大值 500。
aCos = cos(2*pi*50*t);
vCos = aCos.*x1;
sum(vCos)
figure;subplot(3,1,1);plot(x1);subplot(3,1,2);plot(aCos);subplot(3,1,3);plot(vCos);title('check 50Hz');
>> 500
aCos = cos(2*pi*0*t);
vCos = aCos.*x1;
sum(vCos)
figure;subplot(3,1,1);plot(x1);subplot(3,1,2);plot(aCos);subplot(3,1,3);plot(vCos);title('check 0Hz');
>> 1.7808e-13
来做一个全范围测试，写个循环。我们来构建从0Hz到500Hz的501个波（以1Hz为单位逐渐增加）。为什么只到500Hz呢？因为原始信号采样频率是1000Hz（每秒采1000个点），所以我们新构建的波，为了点积求和，必须保持长度一致，一秒也是1000个点；一个周期至少要有一正一负，1000个点最多只能表示到500Hz，再往上就不具有周期性，没意义了。
cs = []; %贮存结果
for i =0:500 %创建不同501个频率的cos波
aCos = cos(2*pi*i*t); %创建当前频率的波
vCos = aCos.*x1; %点积
cs = [cs sum(vCos)]; %求和
end
figure;plot(0:1:500,cs);title('cos 0Hz-500Hz');
cos(2*pi*w1*t).*cos(2*pi*w2*t) = (cos(2*pi*(w1+w2)*t)+cos(2*pi*(w1-w2)*t))/2;第一部分始终是周期波，具有对称性，正负各点互相抵消，点积求和为0。第二部分呢，绝大部分情况下，也是一个周期波，点积求和为0。当且仅当w1==w2时，频率为0，值为1，才可以积累出一个大值。
（b）确定振幅A
正相关显示求和结果为500，为什么是500呢？这个500又怎样对应到原始振幅上呢？
直接看上面推导公式，每个采样点供应了1/2，一共N=1000个采样点，所以最终求和为500。
更一般化，求和得到的值，在最理想（原始波和测试波频率、振幅都一致）的状况下，得到的求和结果就为 A*N/2，A为原始波的振幅，N为采样点。并且推出，采样点频率越高，或者采样时间越长，都会增加样点数从而增高最终的和。
从最终的和反推出原始振幅就很简单了，s=A*N/2 => A=s*2/N;
我们可以试一下，求和值只和样点数与原始振幅相关，和具体的频率无关。
比如30Hz点乘30Hz，点积求和为500：
x1 = cos(2*pi*30*t) ;
aCos = cos(2*pi*30*t);
vCos = aCos.*x1;
sum(vCos)
>> 500
比如300Hz点乘300Hz，点积求和还是500：
x1 = cos(2*pi*300*t) ;
aCos = cos(2*pi*300*t);
vCos = aCos.*x1;
sum(vCos)
>> 500
但是，但是，在两个极端情况，频率为极小值0Hz，或者频率为极大值500Hz时，点积求和却为1000，这又是为什么呢？
x1 = cos(2*pi*0*t) ;
aCos = cos(2*pi*0*t);
vCos = aCos.*x1;
sum(vCos)
>> 1000
x1 = cos(2*pi*500*t) ;
aCos = cos(2*pi*500*t);
vCos = aCos.*x1;
sum(vCos)
>> 1000
这是因为，当极小或极大时，点积之后的波形都变成了水平波，每个样点自己就都是1，不用再凑对子了。多少个样点就是多少和。0Hz的情况下，cos(0)^2 = 1^2 = 1，始终等于1；500Hz的情况下，因为这是当前采样频率下（1000个采样点/秒）能表示的最高频率波，波形就是[1,-1,1,-1....,1,-1]，它的平方也始终是1。在这两个极端情况下，s=A*N => A = s/N;
还有一个小问题，我们的测试波是从0Hz到500Hz，每次增加1Hz。聪明的读者会问，为什么每次增加1Hz呢？可不可以只增加0.1Hz，0.01Hz，这样会测的更准确些吗？我们来试试看。
cs = [];
for i =0:0.1:500 %步长为0.1
aCos = cos(2*pi*i*t);
vCos = aCos.*x1;
cs = [cs sum(vCos)];
end
figure;plot(0:0.1:500,cs);title('cos 0Hz-500Hz');
我们发现，在整数频率上不相关，但是再细分到小数频率上，表现出了各种不可忽略的相关性，这是怎么回事呢？
这是因为，如果测试波的频率不为整数的话，最后总会有个尾巴不是完整周期，其对称性被破坏，无法求和为0。所以测试波必须要有完整的周期。
那么怎么样能把测试精度提升到小数位置呢？我们的例子里，测试的是整一秒的采样，如果我们继续多采样几秒，比如采样2秒钟，那就是2000个采样点，那么所有0.5Hz（例如0.5,1.5,2.5)在一个周期有尾巴，在两个周期正好凑成一个整尾巴，又是一个整周期啦。所以2秒采样时，我们就精确到0.5Hz；3秒采样时，精确到1/3Hz。依此类推，精度=1/采样秒数。
上面的例子，我们想精确到0.1Hz，就要采样10整秒。
Fs = 1000; %采样频率，每秒采样1000点
T = 1/Fs; %每一点就代表1/1000秒
L = 10000; %一共采样10000个点（就是采样10整秒啦）
t = (0:L-1)*T; %每一点对应的时刻（从最开始的0秒，到最后的10秒）
x1 = cos(2*pi*50*t) ; %50Hz信号及采样
figure;plot(t,x1);title('x1 50Hz');
cs = [];
for i =0:0.1:500
aCos = cos(2*pi*i*t);
vCos = aCos.*x1;
cs = [cs sum(vCos)];
end
figure;plot(0:0.1:500,cs);title('cos 0Hz-500Hz 10 seconds');
（c）确定相位b
最后，我们来给原始信号加上相位，然后看看怎样测出这个相位。
Fs = 1000; %采样频率，每秒采样1000点
T = 1/Fs; %每一点就代表1/1000秒
L = 1000; %一共采样1000个点（就是采样一整秒啦）
t = (0:L-1)*T; %每一点对应的时刻（从最开始的0秒，到最后的1秒）
x1 = cos(2*pi*50*t+pi/3) ; %50Hz信号及采样,由于相位关系，整体向左位移pi/3
figure;plot(t,x1);title('x1 50Hz pi/3');
可以看到，由于相位关系，整体向左位移pi/3。
此时再把这个偏移了的信号与50Hz的测试信号点积求和，会得到什么结果呢？
aCos = cos(2*pi*50*t);
vCos = aCos.*x1;
sum(vCos)
figure;subplot(3,1,1);plot(t,x1);subplot(3,1,2);plot(t,aCos);subplot(3,1,3);plot(t,vCos);title('check 50Hz 0pi <-> pi/3');
>> 250.0000
for c =-1:0.01:1
aCos = cos(2*pi*50*t+c*pi);
vCos = aCos.*x1;
cs = [cs,sum(vCos)];
end
figure;plot(-1:0.01:1,cs);title('cos 50Hz -pi:pi');
vCos = aCos.*x1;
sum(vCos)
figure;subplot(3,1,1);plot(t,x1);subplot(3,1,2);plot(t,aCos);subplot(3,1,3);plot(t,vCos);title('check 50Hz 0pi <-> pi/3');
>> 250.0000
可以看到，波形和之前没位移的时候差不多，只是整体下降了1/4。最终求和结果从500减少到了250。
这是因为，cos(2*pi*50*t+b).*cos(2*pi*50*t) = (cos(2*pi*100*t+b)+cos(b))/2;第一部分是周期波，求和之后，各部分互相抵消，和值为0；第二部分是常数，求和之后值为N*cos(b)/2;只有当b=0时，cos(b)取最大值1，和值为500；现在b=pi/3,cos(b)=1/2，和值就只能是250了。
理论上，我们可以找个相位为0的标准波，测一下最大相关和值是多少，以此为参照，通过反余弦，我们就可以得到相位b。但实际上，我们往往无法得到这么个相位为0的标准波（都能控制相位了，还要求什么相位，直接想要什么相位就上什么相位好了）。那怎么办呢？
有两个办法。
方法一：把所有的相位遍历一遍，看看那个相位响应最大。
cs = [];for c =-1:0.01:1
aCos = cos(2*pi*50*t+c*pi);
vCos = aCos.*x1;
cs = [cs,sum(vCos)];
end
figure;plot(-1:0.01:1,cs);title('cos 50Hz -pi:pi');
我们清楚的看到，在0.33pi相位上，表现出了最大的正相关性，所以相位为0.33*pi。
但是遍历的方法，不仅耗时，精度也不高。
方法二：求正弦sin响应值，与余弦cos响应值求比，得出tan，再反推b。
还是这个例子，cos(2*pi*50*t+b).*sin(2*pi*50*t) = (sin(2*pi*100*t+b)+sin(-b))/2;第一部分还是周期波，求和为0；第二部分为常数-sin(b)/2，求和为-N*sin(b)/2;
假设之前那个余弦波测得的响应值为N*cos(b)/2 = v1，现在正弦波测得的响应和值为-N*sin(b)/2=v2;那么有v2/v1 = -sin(b)/cos(b) = -tan(b) => b = atan(-v2/v1)。
并且因为sqrt(cos(b)^2+sin(b)^2) = 1,所以sqrt(v1^2+v2^2)就是0相位时的完全相关求和值。
aCos = cos(2*pi*50*t);
vCos = aCos.*x1;
v1 = sum(vCos)
aSin = sin(2*pi*50*t);
vSin = aSin.*x1;
v2 = sum(vSin)
b = atan(-v2/v1)
b/pi
>> v1 =
250.0000
v2 =
-433.0127
b =
1.0472
b/pi =
0.3333
（d）信号叠加
最后我们再来研究一下，如果原始信号是多（频率）波叠加在一起，还可以用这个相关性方法正确求出各A，w，b吗？比如一个50Hz，一个20Hz的信号加在一起，会不会变成30Hz或者80Hz什么的？相位也互不影响吗？振幅还可以顺利分开来吗？
Fs = 1000; %采样频率，每秒采样1000点
T = 1/Fs; %每一点就代表1/1000秒
L = 1000; %一共采样1000个点（就是采样一整秒啦）
t = (0:L-1)*T; %每一点对应的时刻（从最开始的0秒，到最后的1秒）
x1 = cos(2*pi*50*t) ; %50Hz信号及采样
x2 = 2*cos(2*pi*20*t) ; %20Hz信号及采样，主意它的振幅是2
y = x1+x2;
figure;
subplot(4,1,1);plot(t,x1);title('x1 50Hz');
subplot(4,1,2);plot(t,x2);title('x1 20Hz');
subplot(4,1,3);plot(t,y);title('x1+x2');
cs = [];
for i =0:500
aCos = cos(2*pi*(i)*t);
vCos = aCos.*y;
cs = [cs,sum(vCos)];
end
subplot(4,1,4);plot(0:1:500,cs);title('cos 0Hz-500Hz');
快速傅里叶变换fft，本质上也是求cos sin相关性，只是它把cos放到实数域，sin放到复数域。通过复数的一些对称性加快了计算。将计算复杂度从O(N*N)降低到O(N*log(N))。但是fft必须要输入的样点数为2的次方（2，4，。。。512，1024），不足补零。所以取样时要小心。
（二）调制与解调
有了上面的基础，我们再看调制解调就很简单了。
（a）调制
调制就是把低频波与高频波点乘之后再传输。
Te = 0.01; %每个采样点代表0.01秒
N = 1024; %一共采样1024点，也就是10.24秒
B = 10/(N*Te); %参考频率为10.24秒传输10个周期，也就是0.9766Hz
f1 = B/2; %原始信号为0.4883Hz(当然你可以自己定义）
f0 = 10*B; %载波信号为9.766Hz(当然你可以自己定义）
m = 0.8; %原始信号的振幅为0.8
phi = 0; %载波相位
indice = (0:1:N-1); %各采样点的下标
t = Te*indice; %各采样点的时刻
At = (1+m*cos(2*pi*f1*t+phi)); %原始信号
f = fft(At); %快速傅里叶变换，产生实数部分与虚数部分，分别对应cos,sin
f = abs(f); %求和值
figure;
subplot(2,1,1);plot(indice,At);
subplot(2,1,2);plot(indice,f);
我们可以通过fft观察到，原始信号At，有一个频率为0/10.24 = 0Hz，振幅为1024/1024 = 1的信号；还有一个频率为5/10.24=0.4883Hz，振幅为409.6*2/1024 = 0.8的波。
再来看载波。
C = cos(2*pi*f0*t);
f = fft(C);
f = abs(f);
figure;
subplot(2,1,1);plot(indice,C);title('transport signal');
subplot(2,1,2);plot(indice,f);title('fft of transport signal');
我们看到，载波频率为100/10.24=9.7656Hz。
接着把原始信号与载波信号点乘。
假设之前那个余弦波测得的响应值为N*cos(b)/2 = v1，现在正弦波测得的响应和值为-N*sin(b)/2=v2;那么有v2/v1 = -sin(b)/cos(b) = -tan(b) => b = atan(-v2/v1)。
并且因为sqrt(cos(b)^2+sin(b)^2) = 1,所以sqrt(v1^2+v2^2)就是0相位时的完全相关求和值。
aCos = cos(2*pi*50*t);
vCos = aCos.*x1;
v1 = sum(vCos)
aSin = sin(2*pi*50*t);
vSin = aSin.*x1;
v2 = sum(vSin)
b = atan(-v2/v1)
b/pi
>> v1 =
250.0000
v2 =
-433.0127
b =
1.0472
b/pi =
0.3333
（d）信号叠加
最后我们再来研究一下，如果原始信号是多（频率）波叠加在一起，还可以用这个相关性方法正确求出各A，w，b吗？比如一个50Hz，一个20Hz的信号加在一起，会不会变成30Hz或者80Hz什么的？相位也互不影响吗？振幅还可以顺利分开来吗？
Fs = 1000; %采样频率，每秒采样1000点
T = 1/Fs; %每一点就代表1/1000秒
L = 1000; %一共采样1000个点（就是采样一整秒啦）
t = (0:L-1)*T; %每一点对应的时刻（从最开始的0秒，到最后的1秒）
x1 = cos(2*pi*50*t) ; %50Hz信号及采样
x2 = 2*cos(2*pi*20*t) ; %20Hz信号及采样，主意它的振幅是2
y = x1+x2;
figure;
subplot(4,1,1);plot(t,x1);title('x1 50Hz');
subplot(4,1,2);plot(t,x2);title('x1 20Hz');
subplot(4,1,3);plot(t,y);title('x1+x2');
cs = [];
for i =0:500
aCos = cos(2*pi*(i)*t);
vCos = aCos.*y;
cs = [cs,sum(vCos)];
end
subplot(4,1,4);plot(0:1:500,cs);title('cos 0Hz-500Hz');
可以看到，它们确实是互不影响的。直接套用之前的公式：
（A1*cos(2*pi*w1+b1)+A2*cos(2*pi*w2+b2)).*cos(2*pi*w3+b3) = A1*cos(2*pi*w1+b1).*cos(2*pi*w3+b3)+A2*cos(2*pi*w2+b2).*cos(2*pi*w3+b3);
所以求叠加信号的相关性，和分别求相关性再叠加是一样的，各个信号的计算互不影响。
快速傅里叶变换fft，本质上也是求cos sin相关性，只是它把cos放到实数域，sin放到复数域。通过复数的一些对称性加快了计算。将计算复杂度从O(N*N)降低到O(N*log(N))。但是fft必须要输入的样点数为2的次方（2，4，。。。512，1024），不足补零。所以取样时要小心。
（二）调制与解调
有了上面的基础，我们再看调制解调就很简单了。
（a）调制
调制就是把低频波与高频波点乘之后再传输。
Te = 0.01; %每个采样点代表0.01秒
N = 1024; %一共采样1024点，也就是10.24秒
B = 10/(N*Te); %参考频率为10.24秒传输10个周期，也就是0.9766Hz
f1 = B/2; %原始信号为0.4883Hz(当然你可以自己定义）
f0 = 10*B; %载波信号为9.766Hz(当然你可以自己定义）
m = 0.8; %原始信号的振幅为0.8
phi = 0; %载波相位
indice = (0:1:N-1); %各采样点的下标
t = Te*indice; %各采样点的时刻
At = (1+m*cos(2*pi*f1*t+phi)); %原始信号
f = fft(At); %快速傅里叶变换，产生实数部分与虚数部分，分别对应cos,sin
f = abs(f); %求和值
figure;
subplot(2,1,1);plot(indice,At);
subplot(2,1,2);plot(indice,f);
再来看载波。
C = cos(2*pi*f0*t);
f = fft(C);
f = abs(f);
figure;
subplot(2,1,1);plot(indice,C);title('transport signal');
subplot(2,1,2);plot(indice,f);title('fft of transport signal');
接着把原始信号与载波信号点乘。
Z = At.*C;
figure;
subplot(2,1,1);plot(indice,Z);title('modulated signal');
f = fft(Z);
f = abs(f);
subplot(2,1,2);plot(indice,f);title('fft of modulated signal');
figure;
subplot(2,1,1);plot(indice,Z);title('modulated signal');
f = fft(Z);
f = abs(f);
subplot(2,1,2);plot(indice,f);title('fft of modulated signal');
我们看到，相乘之后波形发生了变换，进一步用fft进行分析，我们发现调制后的信号由3个频率组成，
95/10.24=9.2773Hz，振幅为204.8*2/1024=0.4，相位不明；
105/10.24=10.2539Hz，振幅为204.8*2/1024=0.4，相位不明；
最大的一个：100/10.24=9.7656Hz，振幅为512*2/1024=1，相位不明。
这三个频率都是从哪里来的呢？又有什么关系呢？
At.*C = (1+m*cos(2*pi*f1*t)).*cos(2*pi*f0*t) = cos(2*pi*f0*t)+m.*cos(2*pi*f1*t).*cos(2*pi*f0*t) = cos(2*pi*f0*t)+m*(cos(2*pi*(f0+f1)*t)+cos(2*pi*(f0-f1)*t)/2;
可以从上式发现，包括三个频率f0,f0+f1,f0-f1，以及它们的振幅。
同时可以发现一个规律，除了0Hz信号之外，其它频率f1的信号，载波f0之后，都会分成两个信号f0+f1,f0-f1的叠加。
以上是理想情况，我们再来加噪声看看。
v = 0.08;
N1 = sqrt(v)*randn(1,N)*5;
N2 = sqrt(v)*randn(1,N)*5;
N1 = sqrt(v)*randn(1,N)*5;
N2 = sqrt(v)*randn(1,N)*5;
f = fft(N1);
f = abs(f);
figure;
subplot(2,1,1);plot(indice,N1);title('noise N1');
subplot(2,1,2);plot(indice,f);title('fft noise N1');
f = abs(f);
figure;
subplot(2,1,1);plot(indice,N1);title('noise N1');
subplot(2,1,2);plot(indice,f);title('fft noise N1');
我们通过噪声的fft看到，噪声没有什么特别突出的频率，所以叫噪声。
依葫芦画瓢，对被噪声污染了的信号，我们也看看fft：
Z = At.*C+N1.*cos(2*pi*f0*t)+N2.*sin(2*pi*f0*t);
% Z = At+N1.*cos(2*pi*f0*t)+N2.*sin(2*pi*f0*t);
figure;
subplot(2,1,1);plot(indice,Z);title('noised modulated signal');
% Z = At+N1.*cos(2*pi*f0*t)+N2.*sin(2*pi*f0*t);
figure;
subplot(2,1,1);plot(indice,Z);title('noised modulated signal');
f = fft(Z);
fAbs = abs(f);
subplot(2,1,2);plot(indice,fAbs);
fAbs = abs(f);
subplot(2,1,2);plot(indice,fAbs);
我们看到，即便受到了噪声影响，三个峰值频率位置还是准确的，虽然幅度由于噪声影响发生了变化。
（b）解调
解调就是从高频波里提取原始低频波。
假设我们接收到了受到噪声污染的信号，如上图，我们可以先做下频域上的滤波。因为我知道载波是9.7656Hz（手册，打电话问同事，或者由对称性看出来的。。。），我们就在fft频谱图上，以9.7656Hz为中心（上图中即x=100），左右各保留10个单位，这样既能把原始信号的频率保留下来，由能把噪声的频率尽可能的去除。
filtre = zeros(1, N); %filtre dans le domaine fr¨¦quentiel
filtre(1,100-10:100+10) = 1;
filtre(1,N-100-10:N-100+10) = 1;
fAbsFilter = filtre.*fAbs;
figure;plot(indice,fAbsFilter);
看，是不是清爽多了？
怎样能从调频之后的频率，恢复到原始的低频率信号呢（解调）？再看一下上面的公式：
At.*C = (1+m*cos(2*pi*f1*t)).*cos(2*pi*f0*t) = cos(2*pi*f0*t)+m.*cos(2*pi*f1*t).*cos(2*pi*f0*t) = cos(2*pi*f0*t)+m*(cos(2*pi*(f0+f1)*t)+cos(2*pi*(f0-f1)*t)/2;
我们现在从调频后的高频fft，可以轻松恢复最后一项蓝字部分。把每一部分的频率，减去载波频率f0：
cos(2*pi*(f0-f0)*t)+m*(cos(2*pi*(f0+f1-f0)*t)+cos(2*pi*(f0-f1-f0)*t)/2
= cos(2*pi*0*t)+m*(cos(2*pi*f1*t)+cos(2*pi*f1*t)/2
=1+m*cos(2*pi*f1*t) = At;
filtre = zeros(1, N); %filtre dans le domaine fr¨¦quentiel
filtre(1,100-10:100+10) = 1;
filtre(1,N-100-10:N-100+10) = 1;
fAbsFilter = filtre.*fAbs;
figure;plot(indice,fAbsFilter);
看，是不是清爽多了？
怎样能从调频之后的频率，恢复到原始的低频率信号呢（解调）？再看一下上面的公式：
At.*C = (1+m*cos(2*pi*f1*t)).*cos(2*pi*f0*t) = cos(2*pi*f0*t)+m.*cos(2*pi*f1*t).*cos(2*pi*f0*t) = cos(2*pi*f0*t)+m*(cos(2*pi*(f0+f1)*t)+cos(2*pi*(f0-f1)*t)/2;
我们现在从调频后的高频fft，可以轻松恢复最后一项蓝字部分。把每一部分的频率，减去载波频率f0：
cos(2*pi*(f0-f0)*t)+m*(cos(2*pi*(f0+f1-f0)*t)+cos(2*pi*(f0-f1-f0)*t)/2
= cos(2*pi*0*t)+m*(cos(2*pi*f1*t)+cos(2*pi*f1*t)/2
=1+m*cos(2*pi*f1*t) = At;
以上是phi等于0，原始信号没有相位的情况。给phi赋值pi/3，又会怎样呢？先从理论上看一下相位的影响：
At.*C = (1+m*cos(2*pi*f1*t+b)).*cos(2*pi*f0*t) = cos(2*pi*f0*t)+m.*cos(2*pi*f1*t+b).*cos(2*pi*f0*t) = cos(2*pi*f0*t)+m*(cos(2*pi*(f0+f1)*t+b)+cos(2*pi*(f0-f1)*t-b)/2;
收到蓝字部分的信号后，同样各部分频率减去f0：
cos(2*pi*(f0-f0)*t)+m*(cos(2*pi*(f0+f1-f0)*t+b)+cos(2*pi*(f0-f1-f0)*t-b)/2 = cos(0)+m*(cos(2*pi*f1*t+b)+cos(2*pi*(-f1)*t-b))/2 = cos(0)+m*(cos(2*pi*f1*t+b)+cos(2*pi*f1*t+b))/2 = 1+m*cos(2*pi*ft*t+b)=At;
可以看出，即使有相位，依旧可以通过直接减去载波频率f0，完全恢复。我们在第一部分已经知道，求各个频率的相位是用sin波和cos波本别点积求和，再用atan求出的相位。如果直接用fft呢，其实cos和sin的结果分别保存在实部和虚部，直接用就行了。
下面以有相位的具体例子，结束本教程。本教程力求解释透彻，多有啰嗦，希望对各位看官能有略微帮助。
% modulation
Te = 0.01;
N = 1024;
B = 10/(N*Te);
f1 = B/2;
f0 = 10*B;
m = 0.8;
phi = pi/3;
indice = (0:1:N-1);
t = Te*indice;
At = (1+m*cos(2*pi*f1*t+phi));
f = fft(At);
f = abs(f);
figure;
subplot(2,1,1);plot(indice,At);
subplot(2,1,2);plot(indice,f);
C = cos(2*pi*f0*t);
figure;plot(C);
Z = At.*C;
figure;plot(Z);
f = fft(Z);
f = abs(f);
figure;plot(f);
v = 0.09;
N1 = sqrt(v)*randn(1,N);
N2 = sqrt(v)*randn(1,N);
figure;plot(N1);
f = fft(N1);
f = abs(f);
figure;plot(f);
Z = At.*C+N1.*cos(2*pi*f0*t)+N2.*sin(2*pi*f0*t);
figure;plot(Z);
f = fft(Z);
fAbs = abs(f);
figure;plot(fAbs);
% demodulation
filtre = zeros(1, N); %filtre dans le domaine fréquentiel
filtre(1,100-10:100+10) = 1;
filtre(1,N-100-10:N-100+10) = 1;
figure;plot(filtre);
fFilter = filtre.*fAbs;
figure;plot(fFilter);
ff0 = 101;
ff1 = 96;
ff2 = 106;
rec0 = fAbs(ff0)*cos(2*pi*(ff0-ff0)*100*t/1024-atan2(imag(f(ff0)),real(f(ff0))));
rec1 = fAbs(ff1)*cos(2*pi*(ff0-ff1)*100*t/1024-atan2(imag(f(ff1)),real(f(ff1))));
rec2 = fAbs(ff2)*cos(2*pi*(ff0-ff2)*100*t/1024-atan2(imag(f(ff2)),real(f(ff2))));
figure(1);
subplot(2,1,2);plot((rec0+rec1+rec2)*2/N);title('demodulated signal');
atan2(imag(f(ff0)),real(f(ff0)))/pi %输出3个相位
atan2(imag(f(ff1)),real(f(ff1)))/pi
atan2(imag(f(ff2)),real(f(ff2)))/pi
>> （以pi为单位）
ans =
0.0231
ans =
-0.2749
ans =
0.3055
我们看到，噪声对频率没什么影响，对振幅和相位产生了略微影响，还能接受。
后补：
解调的时候，如果已知原始波的频率，和载波的频率，我们自然可以手动找出调制后的相关频率，然后再手动减去载波频率来解调。
但有的时候，原始波的频率是一段范围，并不固定在某一个点上（干扰什么的，或者像人声，本来就不固定），我们可以直接把收到的信号，重新点乘载波，再做fft。
重新与载波相乘，本质上还是降低信号的频率。我们之前已经看到了，载波可以改变原始信号的频率w1.*wt ->(wt+w1) & (wt-w1)，那么，w1.*wt.*wt -> [(wt+w1) & (wt-w1)].*wt = (2wt+w1) & -w1 & 2wt-w1 & w1。虽然因为对称性，产生了一些原来没有的频率，但原始频率确实能还原回来。比如1Hz的原始波，被10Hz的载波相载，得到调制波9Hz和11Hz，再重新乘以载波，9Hz和11Hz分别还原成1Hz,19Hz,-1Hz,21Hz。（推导公式与之前一样。）
因为我们已经知道波段大概范围，我们只提取出那一小段fft（滤波），再做ifft，x2（因为有一半能量被滤掉了）即可。
dm = Z.*C;
figure;plot(dm);
f = fft(dm);
figure;plot(abs(f));
title('fft after re-multiply tranporter');
filtre = zeros(1,N);
filtre(1,1:10) = 1;
filtre(1,end-10+1:end) = 1;
ff = filtre.*f;
dmm = real(ifft(ff))*2;
figure;plot(dmm);
title('Demodulation first by *T then by ifft');
另外的问题，如果因为传输的延迟，载波的相位也许不为0了。在无法确知载波相位的情况下，又如何恢复原始信号A呢？
原始信号cos(w0*t)，载波cos(w1*t+phi).
调制信号m = cos(w0*t).*cos(w1*t+phi) = 0.5*(cos((w0+w1)*t+phi) + cos((w0-w1)*t-phi));
m.*cos(w1) = 0.5*cos(w0*t)*cos(phi) = 0.5*A*cos(phi);
m.*sin(w1) = 0.5*cos(w0*t)*sin(phi) = 0.5*A*sin(phi);
所以2*sqrt((m.*cos(w1))^2+(m.*sin(w1))^2) = 2*0.5*A = A, 即原始信号。
（全文完 Fou）
At.*C = (1+m*cos(2*pi*f1*t+b)).*cos(2*pi*f0*t) = cos(2*pi*f0*t)+m.*cos(2*pi*f1*t+b).*cos(2*pi*f0*t) = cos(2*pi*f0*t)+m*(cos(2*pi*(f0+f1)*t+b)+cos(2*pi*(f0-f1)*t-b)/2;
收到蓝字部分的信号后，同样各部分频率减去f0：
cos(2*pi*(f0-f0)*t)+m*(cos(2*pi*(f0+f1-f0)*t+b)+cos(2*pi*(f0-f1-f0)*t-b)/2 = cos(0)+m*(cos(2*pi*f1*t+b)+cos(2*pi*(-f1)*t-b))/2 = cos(0)+m*(cos(2*pi*f1*t+b)+cos(2*pi*f1*t+b))/2 = 1+m*cos(2*pi*ft*t+b)=At;
可以看出，即使有相位，依旧可以通过直接减去载波频率f0，完全恢复。我们在第一部分已经知道，求各个频率的相位是用sin波和cos波本别点积求和，再用atan求出的相位。如果直接用fft呢，其实cos和sin的结果分别保存在实部和虚部，直接用就行了。
下面以有相位的具体例子，结束本教程。本教程力求解释透彻，多有啰嗦，希望对各位看官能有略微帮助。
% modulation
Te = 0.01;
N = 1024;
B = 10/(N*Te);
f1 = B/2;
f0 = 10*B;
m = 0.8;
phi = pi/3;
indice = (0:1:N-1);
t = Te*indice;
At = (1+m*cos(2*pi*f1*t+phi));
f = fft(At);
f = abs(f);
figure;
subplot(2,1,1);plot(indice,At);
subplot(2,1,2);plot(indice,f);
C = cos(2*pi*f0*t);
figure;plot(C);
Z = At.*C;
figure;plot(Z);
f = fft(Z);
f = abs(f);
figure;plot(f);
v = 0.09;
N1 = sqrt(v)*randn(1,N);
N2 = sqrt(v)*randn(1,N);
figure;plot(N1);
f = fft(N1);
f = abs(f);
figure;plot(f);
Z = At.*C+N1.*cos(2*pi*f0*t)+N2.*sin(2*pi*f0*t);
figure;plot(Z);
f = fft(Z);
fAbs = abs(f);
figure;plot(fAbs);
% demodulation
filtre = zeros(1, N); %filtre dans le domaine fréquentiel
filtre(1,100-10:100+10) = 1;
filtre(1,N-100-10:N-100+10) = 1;
figure;plot(filtre);
fFilter = filtre.*fAbs;
figure;plot(fFilter);
ff0 = 101;
ff1 = 96;
ff2 = 106;
rec0 = fAbs(ff0)*cos(2*pi*(ff0-ff0)*100*t/1024-atan2(imag(f(ff0)),real(f(ff0))));
rec1 = fAbs(ff1)*cos(2*pi*(ff0-ff1)*100*t/1024-atan2(imag(f(ff1)),real(f(ff1))));
rec2 = fAbs(ff2)*cos(2*pi*(ff0-ff2)*100*t/1024-atan2(imag(f(ff2)),real(f(ff2))));
figure(1);
subplot(2,1,2);plot((rec0+rec1+rec2)*2/N);title('demodulated signal');
atan2(imag(f(ff0)),real(f(ff0)))/pi %输出3个相位
atan2(imag(f(ff1)),real(f(ff1)))/pi
atan2(imag(f(ff2)),real(f(ff2)))/pi
>> （以pi为单位）
ans =
0.0231
ans =
-0.2749
ans =
0.3055
后补：
解调的时候，如果已知原始波的频率，和载波的频率，我们自然可以手动找出调制后的相关频率，然后再手动减去载波频率来解调。
但有的时候，原始波的频率是一段范围，并不固定在某一个点上（干扰什么的，或者像人声，本来就不固定），我们可以直接把收到的信号，重新点乘载波，再做fft。
重新与载波相乘，本质上还是降低信号的频率。我们之前已经看到了，载波可以改变原始信号的频率w1.*wt ->(wt+w1) & (wt-w1)，那么，w1.*wt.*wt -> [(wt+w1) & (wt-w1)].*wt = (2wt+w1) & -w1 & 2wt-w1 & w1。虽然因为对称性，产生了一些原来没有的频率，但原始频率确实能还原回来。比如1Hz的原始波，被10Hz的载波相载，得到调制波9Hz和11Hz，再重新乘以载波，9Hz和11Hz分别还原成1Hz,19Hz,-1Hz,21Hz。（推导公式与之前一样。）
因为我们已经知道波段大概范围，我们只提取出那一小段fft（滤波），再做ifft，x2（因为有一半能量被滤掉了）即可。
dm = Z.*C;
figure;plot(dm);
f = fft(dm);
figure;plot(abs(f));
title('fft after re-multiply tranporter');
filtre = zeros(1,N);
filtre(1,1:10) = 1;
filtre(1,end-10+1:end) = 1;
ff = filtre.*f;
dmm = real(ifft(ff))*2;
figure;plot(dmm);
title('Demodulation first by *T then by ifft');
另外的问题，如果因为传输的延迟，载波的相位也许不为0了。在无法确知载波相位的情况下，又如何恢复原始信号A呢？
原始信号cos(w0*t)，载波cos(w1*t+phi).
调制信号m = cos(w0*t).*cos(w1*t+phi) = 0.5*(cos((w0+w1)*t+phi) + cos((w0-w1)*t-phi));
m.*cos(w1) = 0.5*cos(w0*t)*cos(phi) = 0.5*A*cos(phi);
m.*sin(w1) = 0.5*cos(w0*t)*sin(phi) = 0.5*A*sin(phi);
所以2*sqrt((m.*cos(w1))^2+(m.*sin(w1))^2) = 2*0.5*A = A, 即原始信号。
（全文完 Fou）
