n_samples  = 100; 
f          = 100; % frequency of signal
fs         = 500; % sampling rate
signal     = cos(2*pi*(0:n_samples-1)*f/fs); 

n_samples  = 100; 
f          = 100; % frequency of signal
fs         = 500; % sampling rate
signal     = exp(1i*2*pi*(0:n_samples-1)*f/fs); 

n_samples   = 100; 
qpsk_signal = (floor(2*rand(1,n_samples))-0.5)/0.5 + 1j*(floor(2*rand(1,n_samples))-0.5)/0.5;

nfft        = 4096;
freq_on_bin = 1000 * 1/nfft; 
signal      = exp(1i*2*pi*freq_on_bin*(0:nfft-1));
signal_fft  = fftshift(20*log10(abs(fft(signal,nfft))));
xaxis       = -0.5:1/nfft:0.5-1/nfft; 
figure
plot(xaxis, signal_fft)

% starting with some noise sigma (noise_sigma)
noise_variance = noise_sigma^2; 
noise_signal   = sqrt(noise_variance) .* (randn(1, n_samples));

% starting with some noise sigma (noise_sigma)
noise_variance = noise_sigma^2; 
noise_signal   = sqrt(noise_variance/2) .* (randn(1, n_samples) + 1i*randn(1, n_samples));

% D is your decimation rate
coefs_poly     = reshape(coefs,D,[]);
reg            = zeros(D,size(coefs_poly,2));
v0             = zeros(1,D)';
input_signal   = [1 zeros(1,200)]; % impulse
output_signal  = zeros(1,length(input_signal)/D);

% polyphase decimation filter is "calculate 1 output from D inputs"
out_idx = 1;
for n=1:D:length(input_signal)-D
    v0  = fliplr(input_signal(n:n+D-1)).'; % get D inputs from input signal
    reg = [v0 reg(:,1:size(reg,2)-1)];     % deliver D inputs, shift register values over one

    % calculate 1 output from D inputs
    % apply all polyphase filters to the data
    for k=1:D                        
      output_reg = output_reg + reg(k,:)*coefs_poly(k,:)'; % sum all polyphase paths together
    end
    
    output_signal(out_idx) = output_reg; % set output
    out_idx                = out_idx+1;
end

% coefs are your FIR filter coefficients (firpm, remez, etc)
% I is your interpolation rate

coefs_poly     = reshape(coefs,I,[]);
reg            = zeros(I,size(coefs_poly,2)); % I x filter_length/I
input_signal   = [1 zeros(1,200)]; % impulse
output_signal  = zeros(1,length(input_signal)*I);

% polyphase interpolation filter is "calculate I outputs from 1 input"
out_idx = 1;
for n=1:length(input_signal)
    reg = [input_signal(n) reg(:,1:size(reg,2)-1)]; % deliver 1 input

    % calculate I outputs from the 1 input
    for path=1:I                 
      output_signal(out_idx) = reg * coefs_reshape(path,:)';  % sum all polyphase paths together
      out_idx = out_idx + 1; 
    end
end