take it with a grain of salt but in case someone needs this I think this code works.
clc;
clear all;
close all;
L = 51;
ws1 = 0.4*pi;
wp1 = 0.5*pi;
Ks1 = 1;
Kp = 2;
wo = (ws1+wp1)/2;
G = 1000;
w = [0:(G)]*pi/G;
W = Ks1*(w<=ws1) + Kp*(w>=wp1)
D = (wo>=w);
N=(L-1)/2
m=0:N
k=1:1001
for i = 0:N
F(i+1,:)=sum(arrayfun(@(k) (W(k).^2).*D(k).*cos(w(k).*i), 1:1000))
end
Always happy to be corrected!

