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

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

% QPSK signal matlab code
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;

% FFT matlab code / plot FFT matlab code
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)

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

% complex noise matlab code
% 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));

x = [1 2 3 4 5 6 7]; 
U = 4; % upsample by 4x
x_upsampled = reshape([x; zeros(U-1, length(x))], 1, [])

x = [1 2 3 4 5 6 7]; 
D = 4; % upsample by 4x
x_downsampled = x(1:D:end)

% polyphase decimator matlab code
% 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

% polyphase interpolator matlab code
% coefs are your FIR filter coefficients (firpm, remez, etc)
% I is your interpolation rate
%% INPUT PARAMETERS 
fs = .25; % input sampling rate
fsI = 1;  % output sampling rate
I   = 4;  % interpolation factor
nfft = 10e3;

% cook up some coefficients
[coefs, del] = firpm(143, [0 .025 .1 .5]*2, [1 1 0 0]);

% input signal is two tones, 1 in passband one in stopband
x = exp(1i*2*pi*.025/fs*(0:nfft-1)) + exp(1i*2*pi*.1/fs*(0:nfft-1));

%% MAIN LOOP
coefs_poly     = reshape(coefs,I,[]);
reg            = zeros(1,size(coefs_poly,2)); % I x filter_length/I
input_signal   = x; % 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_poly(path,:)';  % sum all polyphase paths together
      out_idx = out_idx + 1; 
    end
end

%% EXAMPLE PLOTTING
xax  = -.5:1/nfft:.5-1/nfft;
filt_resp  = fftshift(20*log10(abs(fft(coefs/nfft,nfft)))); 
input_resp  = fftshift(20*log10(abs(fft(x/nfft,nfft)))); 
output_resp = fftshift(20*log10(abs(fft(output_signal/nfft,nfft))));
figure
subplot(211)
plot(xax*fsI, filt_resp-max(filt_resp))
hold on
plot(xax*fs, input_resp-max(input_resp))
hold off
legend('filter', 'input')
subplot(212)
hold on
plot(xax*fsI, filt_resp-max(filt_resp))
plot(xax*fsI, output_resp-max(output_resp))
hold off
legend('filter', 'output')

size(input_signal)
size(output_signal)