Hello folks!
Im trying to implement an square root raised cosine FIR pulse shaping filter in verilogA. I have made two versions of this filter, one with coefficients from matlab and one that calculates its own coefficients. Both impulse responses are practically identical. I have also verified the coefficients in matlab by filtering two random I- and Q-signals. They constellation diagram looks like it should, with distinct constellation points only after applying the filter twice, however...

When I use the same coefficients in my verilog model it behaves allmost exactly like a normal raised cosine filter. It's as if the simulator runs the signal twice through the filter. I have tried using my own filter function and the built in numerator/denominator z-transform. Both produce the same result. The input data has been verified and I have tried different settings to no avail.

My questions are as follows, has anyone encountered a problem like this before? Would some kind spirit like to look at my code and try to find the problem?

Many many thanks in advance
//Anders Jakobsson

Here is the code:
//*****************************************************************************
//  ___ ___ ___  ___
// / __| _ \ _ \/ __|
// \__ \   /   / (__
// |___/_|_\_|_\\___|
//
//*****************************************************************************
// FUNCTION: Square root raised cosine filter
// VERSION: 1.0
// AUTHORS: Anders Jakobsson
//
// Description:
//      A square root raised cosine implementation
//
// PARAMETERS:
//      Period        = Symbol period for which the filter is used. [s]
//                default 1u, range [0:inf)
//
//      Gain        = Gain of filter. [N/A]
//                default 1, range (-inf:inf)
//
//      Delay        = Initial delay of filter. [s]
//                default 0, range [0:inf]
//
//*****************************************************************************


`include "constants.h"
`include "discipline.h"


//== Module definition ==
module SRRC2(in,out);
     input in;
     output out;
     voltage in,out;


     //== Parameters ==
     parameter real Period=1u from [0:inf);
     parameter real Gain=1 from (-inf:inf);
     parameter real Delay=0 from [0:inf];
     parameter real Rolloff=0.5 from (0:1];
     parameter integer Oversample=8 from [1:256];
     parameter integer Groupdelay=4 from [1:256];

     
     //== Variables ==
     real T,G,P,OSR,R;
     real tnorm,n1,n2,d,squaresum;
     real taps[1:2*Groupdelay*Oversample+1];
     integer midpoint,undef,i;


     //== Analog section ==
     analog begin
           @(initial_step) begin

                 //Clear variables
                 undef=inf;
                 squaresum=0;
                 G=Groupdelay;
                 P=Period;
                 OSR=Oversample;
                 R=Rolloff;


                 //Set sample period
                 T=Period/OSR;


                 //Define coefficients
                       //Calculate array midpoint index
                       midpoint=G*OSR+1;

                       //set midpiont coefficient ( to h(1) )
                       tnorm=1/OSR;                                    //Calc. tnorm=t/P
                       n1=4*R*cos((1+R)*`M_PI*tnorm);                        //Calc. term one of numerator
                       n2=(1/tnorm)*sin((1-R)*`M_PI*tnorm);                  //Calc. term two of numerator
                       d=`M_PI*(1-pow((4*R*tnorm),2))*sqrt(T);                  //Calc. denominator
                       taps[midpoint]=(n1+n2)/d;                        //Calc. midle impulse response

                       //set other coefficients. only half calculated and mirrored to other half
                       for(i=1;i<=midpoint-1;i=i+1) begin
                             tnorm=(i+1)/OSR;                        //Calc. tnorm=t/P                        
                             n1=4*R*cos((1+R)*`M_PI*tnorm);                  //Calc. term one of numerator
                             n2=(1/tnorm)*sin((1-R)*`M_PI*tnorm);            //Calc. term two of numerator
                             d=`M_PI*(1-pow((4*R*tnorm),2))*sqrt(T);            //Calc. denominator
   
                             //Check for zero denominator
                             if(d!=0) begin
                                   taps[midpoint+i]=(n1+n2)/d;            //Set point and "mirror point"
                                   taps[midpoint-i]=taps[midpoint+i];
                             end else undef=i;
                       end

                       if(undef!=inf) begin                              //If the denominator is zero, the value
                             i=midpoint+undef;                        //is interpolated from the adjacent values
                             taps[i]=( taps[i-1]+taps[i+1] )/2;
                             i=midpoint-undef;
                             taps[i]=( taps[i-1]+taps[i+1] )/2;
                       end


                       //Scale coefficients to some decent value
                       for(i=1;i<=2*G*OSR+1;i=i+1) squaresum=squaresum+pow(taps[i],2);
                       squaresum=5*sqrt(squaresum);
                       for(i=1;i<=2*G*OSR+1;i=i+1) taps[i]=taps[i]/squaresum;

           end

           //Filters the input signal
           V(out) <+ Gain*zi_nd(V(in), taps, {1}, T,T,Delay);
     end
endmodule