用matlab表示周期方波序列（离散时间的）并画出频谱

fs = 100000; %采样频率
T = 1/fs; %采样间隔
N = 1000; %采样点数
t = -0.002:T:0.002;
f = -fs/2 : fs/N : fs/2-fs/N;

% 生成周期方波序列
x = square(2*pi*1000*t,50); 

% 绘制原始信号及其频谱图
figure(1);
subplot(2,1,1);
plot(t,x);
xlabel('Time');
ylabel('Amplitude');
title('Original Signal: Periodic Square Wave');
grid on;
subplot(2,1,2);
sp_x = abs(fftshift(fft(x,N)))/N;
plot(f,sp_x);
xlabel('Frequency');
ylabel('Magnitude');
title('Frequency Spectrum of Original Signal');
grid on;

% 设计滤波器
Wp = 3100*2*pi/(fs*pi); % 通带频率
Wr = 4800*2*pi/(fs*pi); % 阻带频率
Ap = 1; % 通带最大纹波
Ar = 60; % 阻带最小衰减

[N1,wc1] = ellipord( Wp , Wr , Ap , Ar );
[b,a] = ellip(N1,Ap,Ar,wc1,'low'); 

% 绘制滤波器的幅相特性
figure(2);
freqz(b,a,1024,fs);
title('Magnitude and Phase Response of Elliptic Low-pass Filter');

% 使用滤波器进行滤波
y = filter(b,a,x);

% 绘制滤波后的信号及其频谱图
figure(3);
subplot(2,1,1);
plot(t,y);
xlabel('Time');
ylabel('Amplitude');
title('Filtered Signal');
grid on;
subplot(2,1,2);
sp_y = abs(fftshift(fft(y,N)))/N;
plot(f,sp_y);
xlabel('Frequency');
ylabel('Magnitude');
title('Frequency Spectrum of Filtered Signal');
grid on;