正弦信号的频谱分析

最新推荐文章于 2024-08-15 01:55:01 发布

Хайде

最新推荐文章于 2024-08-15 01:55:01 发布

阅读量4.9k

点赞数 4

分类专栏： MATLAB 文章标签： python numpy 开发语言

本文链接：https://blog.csdn.net/weixin_40933496/article/details/126309658

版权

MATLAB 专栏收录该内容

15 篇文章 0 订阅

订阅专栏

在这里插入图片描述

clear,clc,close all;
t_start=tic;
tl=0.1; 
f1=20; 
f2=20.02;
fs=60;%抽样频率
N=60;%抽样点数




O0=0;%函数的起始相位角弧度制
cnt1=1;
fs=50; 
for N=10:30:1000 
    disp(strcat(strcat('第',num2str(cnt1)),'次运算'));%Amp_max1从12列之后的是找到的极大值的索引值,inf是出错补的值
   Amp_max1(cnt1,:)=spect_function(f1,fs,N,O0,tl);
    cnt1=cnt1+1;
end
 close all;
N=200;
cnt2=1
for fs=50:50:2500
    disp(strcat(strcat('第',num2str(cnt2)),'次运算'));
    Amp_max2(cnt2,:)=spect_function(f1,fs,N,O0,tl);
    cnt2=cnt2+1;
end
close all;
Amp_max1=Copy_of_spect_function(f1,f2,fs,N,O0,tl);
cnt1=1;
fs=50;
 for N=500:500:5000%N=10:30:1000%
    disp(strcat(strcat('第',num2str(cnt1)),'次运算'));
     Amp_max3(cnt1,:)=Copy_of_spect_function(f1,f2,fs,N,O0,tl);
    cnt1=cnt1+1;
end

close all
N=60000;
cnt2=1
for fs=50:50:1000
    disp(strcat(strcat('第',num2str(cnt2)),'次运算'));
    Amp_max4(cnt2,:)=Copy_of_spect_function(f1,f2,fs,N,O0,tl);
    cnt2=cnt2+1;
end
 close all;
 t_end=toc(t_start);
 disp(strcat('程序的运行时间是:',num2str(t_end)));
被调函数1:
function [ Amp_max ] = spect_function( f1,fs,N ,O0,t1)
n=0:N-1;
x=cos(2*pi*f1*n/fs+O0);
figure;
t=0:0.0001:(N-1)/fs;
subplot(311);
plot(t,cos(2*pi*f1*t+O0));title('cos(2*pi*f1*t+O0)');grid on;
subplot(312)
%f = (0:N-1)*fs/N;
% plot(n*fs/N,x,'*');hold on;plot(n*fs/N,x)
plot(n/fs,x,'*');hold on;plot(n/fs,x);xlabel('t/s');
title(strcat(strcat(strcat('f=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N))));
grid on;
X = fftshift(fft(x./N)); 
f = (-N/2:N/2-1)*(fs/N);
subplot(313)
plot(f,abs(X)); 
hold  on;X_f=abs(X);
X_f_max_index=find(X_f==max(X_f));
plot(f(X_f_max_index),X_f(X_f_max_index),'*');
xlabel('f/Hz');ylabel('幅度')
title(strcat(strcat(strcat('f=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N))));
grid on;
hold off;
Amp_max=[f(X_f_max_index),X_f(X_f_max_index)];
Amp_max=[f1,fs,N,Amp_max];
[M1,N1]=size(Amp_max);
Amp_max(1,8)=0;
while N1<7
    Amp_max=[Amp_max,inf];
    [M1,N1]=size(Amp_max);
    Amp_max(1,8)=1;
end
if N1>7
    Amp_max=Amp_max(1:7); 
    Amp_max(1,8)=-1;
end
disp(strcat(strcat(strcat('f=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N))));
disp(strcat(strcat('幅度最大值对应的频率是:',num2str(f(X_f_max_index))),strcat('   幅度最大值是:',num2str(X_f(X_f_max_index)))));
s1=sum(x.^2);
s2=sum(abs(fft(x)).^2)/N;
Amp_max(1,9)=s1; 
Amp_max(1,10)=s2; 
flag=0;
if abs(s1-s2)<=10^-5
    flag=1;
    disp(strcat(strcat(strcat('时域能量:',num2str(s1)),strcat(',频域能量:',num2str(s2))),'满足帕塞瓦尔定理'));
else
    disp(strcat(strcat(strcat('时域能量:',num2str(s1)),strcat(',频域能量:',num2str(s2))),'不满足帕塞瓦尔定理'));
end
Amp_max(1,11)=flag;
disp('----------------------------------');

if(length(X_f_max_index)~=2)
    if length(X_f_max_index)>2
        X_f_max_index=X_f_max_index(1:2);
    end
    while length(X_f_max_index)<2
        X_f_max_index=[X_f_max_index,inf];
    end
end
Amp_max=[Amp_max,X_f_max_index];%加上最大值的索引k值
% pause(tl);
pause(t1);
end
被调函数2：
function [ Amp_max ] = Copy_of_spect_function( f1,f2,fs,N,O0,t1 )
n=0:N-1;
x=cos(2*pi*f1*n/fs+O0)+cos(2*pi*f2*n/fs+O0);
figure;
t=0:0.0001:(N-1)/fs;
subplot(311);
plot(t,cos(2*pi*f1*t+O0)+cos(2*pi*f2*t+O0));title('原始信号cos(2*pi*f1*t+O0)+cos(2*pi*f2*t+O0)');grid on;
subplot(312)
plot(n/fs,x,'*');hold on;plot(n/fs,x)
xlabel('t/s');
title(strcat(strcat('f2=',num2str(f2)),strcat(strcat(strcat(' f1=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N)))));
%axis([-inf inf -1 1]);
grid on;
X = fftshift(fft(x./N)); 
f = (-N/2:N/2-1)*(fs/N);
subplot(313)
plot(f,abs(X)); 
hold  on;
X_f=abs(X);
X_f_max_index=find(X_f==max(X_f));
plot(f(X_f_max_index),X_f(X_f_max_index),'*');
xlabel('f/Hz')
ylabel('幅度')
%axis([19 ,19.4,-inf inf]);
title(strcat(strcat('f2=',num2str(f2)),strcat(strcat(strcat(' f1=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N)))));
grid on;hold off;figure;subplot(311)
plot(f,abs(X));%画双侧频谱幅度图
hold  on;X_f=abs(X);
X_f_max_index=find(X_f==max(X_f));
plot(f(X_f_max_index),X_f(X_f_max_index),'*');
xlabel('f/Hz')
ylabel('幅度')
title(strcat(strcat('f2=',num2str(f2)),strcat(strcat(strcat(' f1=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N)))));
grid on;
hold off;
subplot(312)
plot(f,abs(X)); 
hold  on;
X_f=abs(X);
X_f_max_index=find(X_f==max(X_f));
plot(f(X_f_max_index),X_f(X_f_max_index),'*');
xlabel('f/Hz')
ylabel('幅度');axis([-20.1 ,-19.9,-inf inf]);
title('负频谱');grid on;hold off;subplot(313)
plot(f,abs(X));%画双侧频谱幅度图
hold  on;X_f=abs(X);
X_f_max_index=find(X_f==max(X_f));
plot(f(X_f_max_index),X_f(X_f_max_index),'*');
xlabel('f/Hz');ylabel('幅度')
axis([19.9 ,20.1,-inf inf]);
title('正频谱');grid on;hold off;
Amp_max=[f(X_f_max_index),X_f(X_f_max_index)
];
Amp_max=[f1,fs,N,Amp_max];
[M1,N1]=size(Amp_max);
Amp_max(1,8)=0;
while N1<7
    Amp_max=[Amp_max,inf];%异常处理机制
    [M1,N1]=size(Amp_max);
    Amp_max(1,8)=1;
end
if N1>7
    Amp_max=Amp_max(1:7);%异常处理机制
    Amp_max(1,8)=-1;
end
disp((strcat(strcat('f2=',num2str(f2)),strcat(strcat(strcat(' f1=',num2str(f1)),strcat(',fs=',num2str(fs))),strcat(',N=',num2str(N))))));
disp(strcat(strcat('幅度最大值对应的频率是:',num2str(f(X_f_max_index))),strcat('   幅度最大值是:',num2str(X_f(X_f_max_index)))));
s1=sum(x.^2);
s2=sum(abs(fft(x)).^2)/N;
Amp_max(1,9)=s1;%信号时域能量
Amp_max(1,10)=s2;%信号频域能量
flag=0;
if abs(s1-s2)<=10^-5
    flag=1;
    disp(strcat(strcat(strcat('时域能量:',num2str(s1)),strcat(',频域能量:',num2str(s2))),'满足帕塞瓦尔定理'));
else
    disp(strcat(strcat(strcat('时域能量:',num2str(s1)),strcat(',频域能量:',num2str(s2))),'不满足帕塞瓦尔定理'));
end
Amp_max(1,11)=flag;
disp('----------------------------------');
pause(t1);
end