spectrumanalyse.m

频谱分析

%2006.3.5
%1、产生40KHZ,70KHZ的信号，分析其频谱分布，并分析原因。
%2、从画出的功率谱中，我们得出40KHZ、70KHZ交替发射时，功率谱大约在55KHZ处出现最大值
  %程序中通过FFT（得出X[K]）和傅立叶级数（得出C）比较，解释频率出现上述分布的原因
  %并给出任意两种频率组合时傅立叶级数的通式（见fourier.m）
 %3、出现的问题：FFT时采样点数不一样，功率谱出现很大差别（目前尚未解决）
clear all;
fs=5000;%采样频率
%a1：只有70KHZ的高电平段
%a2：只有40KHZ的高电平段
%a：叠加后的信号

N=50;%50khz
N1=36;%70khz
N2=63;%40khz
a=zeros(2500,1);
k=1;
j=1;
T0=11/280;%周期

for k=1:12
    for j=1:N1
        a(36+198*(k-1)+j)=1;
    end
end
a1=a;

for k=1:12
    for j=1:N2
        a(135+198*(k-1)+j)=1;
    end
end
a2=a;

for k=1:12
    for j=1:N1
        a2(36+198*(k-1)+j)=0;
    end
end

%用fft画各次谐波的幅度（c=b.*conj(b)/2500画出的才是功率谱）
b=fft(a);
c=b.*conj(b);
c=sqrt(c);

%找峰值点（粗略认为是谐波所在位置），目的是和傅立叶级数求出的谐波幅度做比较，
%   看同次谐波的幅度是否相等
i=1;
for k=2:2499
    if(c(k)>c(k-1)&&c(k)>c(k+1))
        cc(i)=c(k);
        i=i+1;
    else
    end
end    
  
f=5000000*(1:2500)/2500;
f1=5000000*(1:1000)/1000;

b1=fft(a1);
a1=a1(1:1000);
b1t=fft(a1,1000);
c1=b1.*conj(b1);
c1=sqrt(c1);
c1t=b1t.*conj(b1t);%采样点为5个周期时，比较与采全屏幕时得出的频率有什么不同
c1t=sqrt(c1t);

i=1;
for k=2:999
    if(c(k)>c(k-1)&&c(k)>c(k+1))
        cc1(i)=c1(k);
        i=i+1;
    else
    end
end  
i=1;
for k=2:999
    if(c(k)>c(k-1)&&c(k)>c(k+1))
        cc1t(i)=c1t(k);
        i=i+1;
    else
    end
end  

for k=1:86
    ccc(k)=cc1(k)/cc1t(k);
end
    

b2=fft(a2);
c2=b2.*conj(b2);
c2=sqrt(c2);

figure(1)
subplot(321)
plot(a1);
axis([1 2500 -2 2]);
grid on;
subplot(322)
plot(f,c1);
title('2500点快速傅立叶变换');
hold on;
plot(f1,c1t,'r');
axis([0 260000 0 600]);
grid on;
subplot(323)
plot(a2)
axis([1 2500 -2 2]);
grid on;
subplot(324)
plot(f,c2);
axis([0 260000 0 600]);
grid on;
subplot(325)
plot(a);
axis([1 2500 -2 2]);
grid on;
subplot(326)
plot(f,c);
axis([0 260000 0 600]);
grid on;


 %两种不同频率组合时各次谐波的幅度表达式，其中
 %      a1、a2、b1、b2分别表示高电平的起始和终止位置
 % a(n)=(sin(2*pi*n*a2)-sin(2*pi*n*a1)+sin(2*pi*n*b2)-sin(2*pi*n*b1))/(n*pi);
 % b(n)=(cos(2*pi*n*a1)-cos(2*pi*n*a2)+cos(2*pi*n*b1)-cos(2*pi*n*b2))/(n*pi);
 % c=sqrt(a.^2+b.^2)
clear all;
for n=1:10
    a(n)=(sin(2*pi*n*72/198)-sin(2*pi*n*36/198)-sin(2*pi*n*135/198))/(n*pi);
    b(n)=(cos(2*pi*n*36/198)-cos(2*pi*n*72/198)-1+cos(2*pi*n*135/198))/(n*pi);
end
c=sqrt(a.^2+b.^2);

for n=1:10
    c1(n)=sqrt(2-2*cos(36*2*n*pi/198))/(n*pi);
    c2(n)=sqrt(2-2*cos(63*2*n*pi/198))/(n*pi);
end
figure(10)
subplot(311)
plot(10*c1,'diamond');
title('傅立叶级数');
grid on;
subplot(312)
plot(10*c2,'diamond');
grid on;
subplot(313)
plot(1000*c);
grid on;


if(0)
figure(2)
subplot(311)
plot(f1,c11,'diamond');
grid on;
%axis([0 100000 0 500]);
subplot(312)
plot(c12);
grid on;
%axis([0 100000 0 500]);
subplot(313)
plot(c13);
grid on;
%axis([0 100000 0 500]);
%end
figure(3)
plot(c111,'diamond');
hold on;
plot(c112,'r');
hold on;
plot(c113,'b');
grid on;

%对占空比为1/2的信号进行分析
for k=1:2:48
    for j=1:N
        a(50*k+j)=1;
    end
end
b=fft(a);
c=b.*conj(b)/2500;
f=5000000*(1:2500)/2500;
figure(1)
subplot(211)
plot(a);
axis([1 2500 -2 2]);
grid on;
subplot(212)
plot(f,c,'diamond');
axis([0 300000 0 300]);
grid on;
end