How hard is this combinatorial optimisation problem?
This problem can be solved with dynamic programming in pseudo-polynomial time (proof below). Therefore, it is not possible to show that this problem is strongly NP-hard (unless P=NP).

First, let's restate the problem:
Given: values $N$ and $T$ and positive integer intervals $R_1$, $R_2$, $\ldots$, and $R_n$
Output: the largest possible value of $\sum_{i=1}^nf(x_i)$ (where $f$ is defined in terms of $N$ as described above) such that the variables $x_1, x_2, \ldots, x_n$ satisfy all of the following:

$x_i \in R_i$ for $i = 1,2,\ldots,n$
$\sum_{i=1}^nx_i \le T$
for all pairs $i,j$, if $\text{median}(R_i) \le \text{median}(R_j)$ then $\frac{\text{median}(R_i)}{x_i} \le \frac{\text{median}(R_j)}{x_j}$

We can assume without loss of generality that the $R_i$s are sorted by median. That is, assume that $\text{median}(R_1) \le \text{median}(R_2) \le \ldots \le \text{median}(R_n)$. If we do, the final condition can be simplified. It is both necessary and sufficient that $\frac{\text{median}(R_1)}{x_1} \le \frac{\text{median}(R_2)}{x_2} \le \ldots \le \frac{\text{median}(R_n)}{x_n}$.
(Note: this is not precisely correct. Assuming no intervals have the same median, this is correct, but if multiple intervals have the same median it is actually also necessary that the $x_i$ values chosen for those intervals have to all be the same. This makes things more complicated without really making the problem any harder, so I'm going to ignore this possibility for the rest of this answer.)

With that done, we can define the subproblems that we will use for our dynamic programming solution.
Given: all the inputs for the problem as well as (1) an index $k$ with $1 \le k \le n$, (2) a threshold $t$ with $0 \le t \le T$, and (3) a value $X \in R_k$
Output: the largest possible value of $\sum_{i=1}^kf(x_i)$ such that the variables $x_1, x_2, \ldots, x_k$ satisfy all of the following:

$x_i \in R_i$ for $i = 1,2, \ldots, k$
$x_k = X$
$\sum_{i = 1}^kx_i \le t$
for $i = 1,2, \ldots, k-1$, it is the case that $\frac{\text{median}(R_i)}{x_i} \le \frac{\text{median}(R_{i+1})}{x_{i+1}}$

Let $S[k][t][X]$ denote the output of this subproblem.

Note that the overall problem's answer can be expressed in terms of several subproblems: the overall answer is $\max_{X \in R_n}S[n][T][X]$.
Note also that the number of subproblems is small (polynomial in the size of the numbers in the input of the problem). In particular, $k$ takes on $n$ possible values, $t$ takes on $T+1$ possible values, and $X$ takes on $\sum_{i=1}^n|R_i|$ possible values. Each of these is polynomial in the numbers used as inputs to the problem. Therefore, the overall number of subproblems is also polynomial in the numeric value of the input. Thus, provided we can show that each subproblem can be solved in pseudo-polynomial time, it will also be the case that all the subproblems together can be solved in pseudo-polynomial time, and therefore that the overall problem can be solved in pseudo-polynomial time.

In order to actually solve the subproblems, we will use dynamic programming. In other words, we will build up a table of values for $S[k][t][X]$. We can compute these subproblems in order primarily by $k$. In other words, we can assume that by the time we are computing $S[k][t][X]$, the value of $S[k'][t'][X']$ has already been computed provided $k' &lt; k$.
So how do we actually compute these values?
The base case, where $k=1$, is easy. $S[1][t][X]$ is defined as the largest possible value of $\sum_{i=1}^kf(x_i) = f(x_1)$ where $x_1 = X \le t$. Therefore, provided $X \le t$, we have that $S[1][t][X] = f(X)$. Otherwise, no such value is possible and so we can set $S[1][t][X]$ to $-\infty$ to indicate that no real value is achievable.
Otherwise, let's say we're trying to compute $S[k][t][X]$ for some $k &gt; 1$. This value is defined as the largest possible value of $\sum_{i=1}^kf(x_i)$ subject to the following constraints:

$x_i \in R_i$ for $i = 1,2, \ldots, k$
$x_k = X$
$\sum_{i = 1}^kx_i \le t$
for $i = 1,2, \ldots, k-1$, it is the case that $\frac{\text{median}(R_i)}{x_i} \le \frac{\text{median}(R_{i+1})}{x_{i+1}}$

Since we know $x_k = X$, the value we are maximizing is $\sum_{i=1}^kf(x_i) = f(x_k) + \sum_{i=1}^{k-1}f(x_i) = f(X) + \sum_{i=1}^{k-1}f(x_i)$. This is the same as just maximizing $\sum_{i=1}^{k-1}f(x_i)$. What conditions do the variables $x_1, \ldots, x_{k-1}$ have to satisfy? We can reframe the above conditions as follows:

$x_i \in R_i$ for $i = 1,2, \ldots, k-1$
$\sum_{i = 1}^{k-1}x_i \le t - X$
for $i = 1,2, \ldots, k-2$, it is the case that $\frac{\text{median}(R_i)}{x_i} \le \frac{\text{median}(R_{i+1})}{x_{i+1}}$
it is the case that $\frac{\text{median}(R_{k-1})}{x_{k-1}} \le \frac{\text{median}(R_{k})}{X}$

By the final condition, only some of the possible values of $x_{k-1}$ are actually permissible values. In particular, we have the condition that $\frac{\text{median}(R_{k-1})}{x_{k-1}} \le \frac{\text{median}(R_{k})}{X}$, which can be rearranged to be $x_{k-1} \ge X \times \frac{\text{median}(R_{k-1})}{\text{median}(R_{k})}$. Thus, $x_{k-1}$ must be some value from $R_{k-1} \cap [X \times \frac{\text{median}(R_{k-1})}{\text{median}(R_{k})}, \infty)$. Call this other interval $R_{k-1}'(X)$
Suppose we fix some value $Y \in R_{k-1}'(X)$ and add the constraint that $x_{k-1} = Y$. At that point, we're trying to maximize $\sum_{i=1}^{k-1}f(x_i)$ subject to the following:

$x_i \in R_i$ for $i = 1,2, \ldots, k-1$
$x_{k-1} = Y$
$\sum_{i = 1}^{k-1}x_i \le t - X$
for $i = 1,2, \ldots, k-2$, it is the case that $\frac{\text{median}(R_i)}{x_i} \le \frac{\text{median}(R_{i+1})}{x_{i+1}}$

But this is just the definition of $S[k-1][t-X][Y]$.
It should be pretty clear at this point that $S[k][t][X]$ is just the maximum over $Y \in R_{k-1}'(X)$ of $f(X) + S[k-1][t-X][Y]$. Note: if $R_{k-1}'(X)$ is empty then no choice of values is possible and so we again set $S[k][t][X]$ to $-\infty$. Similarly, we can assume that if $t-X &lt; 0$ then the value $S[k-1][t-X][Y]$ is just defined to be $-\infty$ by default. (Note that this formula appropriate handles the case where some of the $S[k-1][t-X][Y]$ values are $-\infty$).
Anyway, assuming you deal with all the edge cases, this shows that you can compute $S[k][t][X]$ in terms of the previous values $S[k'][t'][X']$ with $k' &lt; k$. The runtime of this computation is linear in the size of $R_{k-1}'(X)$, which is at most $|R_{k-1}|$. Therefore, the computation of each $S[k][t][X]$ is pseudo-polynomial in the input, implying (as previously argued) that the computation of all the subproblems together also runs in pseudo-polynomial time.

EDIT: I was asked to address two specific issues in the comments. First of all, what happens if the function $f$ can depend on which interval $R_i$ it is taken from. And second, what happens if not every interval has to contribute a value (i.e. if you are allowed to skip an $x_i$). I believe the problem is still solvable with the same pseudo-polynomial time dynamic programming approach. I also mentioned earlier that the solution above ignores the possibility that multiple ranges share the same medians. I have not edited the solution above, but to address these concerns, I'm adding a coded up version of my algorithm (using python), fixed to address all these issues.
I assume that the input consists of a value T (representing $T$) and a list Rs, whose elements represent the ranges $R_i$ and which can be looped over. I assume that the elements of Rs are in sorted order by median and that there exists a function median that can identify the median of any such range. I also assume that there exists a function f which takes as input two values: an index i and a value x. f(i, x) is used to represent $f_i(x)$.
def solve(T, Rs):
  if len(Rs) == 0:
    return 0

  # S is going to be the subproblem lookup table. S[k][t][X] will denote
  # the largest possible output value using xs only from ranges Rs[0]
  # through Rs[k] (with each range contributing at most one x), such that
  # the sum of the xs is at most t, and such that range Rs[k] does contribute
  # a value, and the value it contributes is X.
  S = []

  # base case
  layer = [{X: float(-inf) for X in Rs[0]} for t in range(T+1)] # default to -inf
  S.append(layer)
  for t in range(T+1):
    for X in Rs[0]:
      if X &lt;= t:
        S[0][t][X] = f(0, X)

  # inductive case
  for k in range(1, len(Rs)):
    layer = [{X: float(-inf) for X in Rs[0]} for t in range(T+1)] # default to -inf
    S.append(layer)
    for t in range(T+1):
      for X in Rs[k]:
        # the minimum necessary condition is that X &lt;= t
        if X &lt;= t:
          # if no other range contributes an x, then the smallest possible
          # value is f(k, X)
          layer[t][X] = f(k, X)
          # if some other range does contribute an x, then there is a largest
          # index range which does so; the index of that range is some i &lt; k
          # and we can loop over the possibilities
          for i in range(k):
            # we can also loop over the possible values of what x that range
            # contributed
            for Y in Rs[i]:
              # in order for this to be valid, the following must hold:
              if median(Rs[i]) / Y &lt;= median(Rs[k]) / X:
                # also, if the medians are the same, we also have the
                # constraint that X=Y; therefore, this case is only valid
                # if the following holds
                if median(Rs[i]) != median(Rs[k]) or X == Y:
                  # the best possible value in this case is the following
                  value = f(k, X) + S[i][t-X][Y]
                  # if this is better than anything so far, use this value
                  S[k][t][X] = max(S[k][t][X], value)

  # Now lets actually get the overall answer. One possibility is to use no xs
  # from any range:
  best = 0

  # otherwise, there will be some largest index range which contributes an x
  # and it will contribute some specific X
  for k in range(len(Rs)):
    for X in Rs[k]:
      # the largest possible value in that case subject to all the constraints:
      value = S[k][T][X]
      # we can use that value instead of the best one saved if it's better
      best = max(best, value)

  return best