What should be done if the mask-values are too small? According to the research paper Multidirectional Scratch Detection and Restoration in Digitized Old Images, we have,

$$H(u, v) = \frac{1}{1 + 0.414 {. \sqrt[{2n}]{\frac {u^*}{D_u}+\frac {v^*}{D_v}}}} \tag{1}$$

where:

\begin{align}
u^* &= \cos \theta . (u + t_x) + \sin \theta . (v + t_y)\\
v^* &= -\sin\theta . (u + t_x) + \cos \theta . (v + t_y)\\
t_x &= \mbox{center}_x \times \cos \theta\\
t_y &= \mbox{center}_y \times \sin \theta
\end{align}

N.B. I have used an 128x128 8-bit-indexed gray-scale image and a 32x32 mask. I padded the mask. Then converted both the image and the mask to complex 2d array. Then I applied Fourier transform to both of them. Then I multiplied them. Then did the Inverse transform.

I have tested with values:

$\theta=0.9$ and $radian=0.9$
$D_u=2, D_v=2$
$CenterX=16, CenterY=16$ and $CenterX=-16, CenterY=-16$

The results are not coming as expected.

The convolution operation is generating a complete black image. My guess is, the kernel values are too small.

Is there any technique to magnify these data?

P.S I have actually corrected formula $(1)$

A few of the stuff that appears in this response were initially packed together in a compact comment. Inevitably, in this answer, I would have to provide code, which I wouldn't like in this case. You have to get what is going on. So, I am going to try and talk a little bit more "around" the code and hopefully by the end of this, the code will have written itself.

In just a few words: $H$ is your frequency response in the frequency domain, it is based on a bandpass pattern but it is supposed to be applied in a directional way. So, build a bandpass and rotate it.

Now then...Bandpass "pattern". You need something that goes up and then down and when it goes up it stays up as flat as possible. Take a plot of some $y_1 = \frac{x}{a}$. Well, obviously, that's a line with a slope of $\frac{1}{a}$ and it crosses zero at zero. Nothing fancy going on there.

Now, offset that zero a bit to the right, so $y_2 = \frac{b-x}{a}$. What did you just do? You brought zero bang on on $b$ but also you introduced some NEGATIVE values for $x>b$.

To see why this will come handy, do a $y_3=(\frac{b-x}{a})^2$. This is a parabola. It goes down and then up. It is also curved, let's remove this curve by undoing the square with a root, so $y_4=\sqrt{(\frac{b-x}{a})^2}$. What's going on now? We have two lines, one that arrives down to $b$ and one that departs upwards from $b$. Let's now reverse these trends.
$y_5=\frac{1}{\sqrt{(\frac{b-x}{a})^2}}$.

What is any $u=\frac{1}{v}$? A hyperbola. Remember those two lines that meet at $b$? If you were to substitute $v$ for $y_5$, then, when $y_5$ goes down, towards $b$, it tends to smaller values, therefore, $u$ goes upwards!. And when $y_5$ departs from $b$ towards increasing values then $y_5$ goes downwards.

What happens bang on on $b$? That's when $y_5$ hits zero and if $v$ hits zero then $u$ goes to infinity. We don't like that. So, let's make that $v$ a sum, to control the denominator and ensure it never hits zero. $y_6 = \frac{1}{1+\sqrt{(\frac{b-x}{a})^2}}$. So now, when that root term goes to zero, $y_6$ stays at $1$.

Let's deal with that stay flat business. Remember that $\sqrt{q^2}$ trick we did to get rid of possible negative values of $q$? Well, we could have done this with $|q|$, which means taking the absolute value of some $q$. But this power raising has an added benefit. So, in general, if you raise to some power, you increase the slope of a function. That is, how sharply the curve evolves. Take some $q^2$ and $q^4$, for any two $q, q+\delta_q$ where $\delta_q > 0$, then, $q^4$ grows faster than $q^2$.

BUT!, raising to a power has another property as well, raising large numbers to a power makes them larger and raising small numbers ($ < 1$) to a power makes them...even smaller!. Why is this important for us? Because as our $y_2$ approaches its $b$ limit, it grows smaller and smaller and smaller and we can accelerate this decrease by raising to some power. Therefore, we now go $\sqrt{q^n}$ and use $n$ to control how fast that function will approach small values and stay small leading to that flatness characteristic.

Now, here is the thing with multiplication, if you raise some $q$ to an odd power, you retain the original sign of $q$ and if that was negative then $\sqrt(q)$ will become complex. We don't like that in this case. So, what do we do? We ensure that some $q$ is raised to an even power by multiplying by two. All together now:

$$y_7 = \frac{1}{1+\sqrt{(\frac{b-x}{a})^{2n}}}$$

OK, fine, how do I apply $y_7$ as a filter to a signal? Pretend that your $x$ is now frequency, mirror $y_7$ around some $x_{\frac{N}{2}}$ (where $N$ is the full length of your impulse response), take its inverse Fourier Transform, throw away the complex part (which should have been zero anyway) and treat what you get back as your impulse response. That is your $h$.

Hold on, if $x$ now represents frequency then what are the $b,a,n$? Consider how $y_7$ goes up towards something and then comes down from that something. So, $b$, is your Center Frequency. The centre frequency of the bandpass filter. $a$ is now your bandwidth, that is how fat that lobe of values that go up and down is. And finally, $n$ is now your filter's degree, that is, how flat it tends to be. In fact, if you raise $n$ to some absurd value (which, remember, it will be doubled), then it will look very much like a square pulse which would be an ideal filter. But, the more your crank up $n$ the more Fourier coefficients you need to synthesize it accurately, and then the impulse response gets longer (that $N$ from earlier) and...anyway, just keep that in mind, it's not important for now.

WHAT THE?!?!? HOW DOES THIS APPLY TO MY ORIGINAL QUESTION?!?!!?!?!?!

We are very close now, all that we have to do is bring in another dimension, so that our $y_7$ now gets a $y$ as well. (I hope using $y$ and $y_7$ is not confusing, they are just symbols)

I am not going to go all over this again because we need to focus on what is really important here. When you bring in another dimension you can now define direction. I am not going to go through the full two dimensional Fourier spectrum explanation but I will just say this. When you look at the frequency spectrum of an image, the low frequencies are at the centre and the high frequencies are at the periphery (by convention).

Therefore, filters that treat all directions equally look like cakes or donuts. It's not hard to imagine why. Take our bandpass profile $y_7$ and mentally rotate it about the end where its zero lies. This will give you a ring (or doughnut) about the centre of the image. Pick any direction, what you see is the bandpass filter and its symmetric image, remember earlier that I said, "...mirror it around some $x_{\frac{N}{2}}$..."?.

But, hey, we now have the luxury to filter along a specific direction and not treat all of them equally. So, let's take our filter profile and apply it to just a specific direction.

How are we going to do this?

Remember our $y_7$? This is effectively an $f(x)$, a function of some $x$ variable which we can assign any meaning we like. It could be some generic $x$, it could be a frequency axis segment. But, at the end of the day, what we did there was to scan $x$ and assign a $y$ value, so the height of $y$ depends on the value of $x$.

What we are going to do now, is scan the frequency plane and assign a height to it, so imagine your $H$ as being a height map. But! We want to be able to control for the direction which translates to rotation.

So, while I am scanning some $u,v$ point on the frequency plane, I am assigning its height on some $u^*, v^*$ which are the rotated coordinates, around the centre $t_x, t_y$.

At this point, I have to say, I can see why the authors want to call $D_v, D_h$ their cutoff frequencies and $tx,ty$ their "sub-band centre frequency" but from the viewpoint of the above explanation, $D_v, D_h$ would correspond to the bandwidth and $tx,ty$ would correspond to the cutoff frequency.

Now, you work in C# and you probably have very good reasons for this but my quick test was performed in Octave (which is very similar to MATLAB). It is a piece of FOSS and that's where you could have just done a imagesc(H) and looking at your results we could more or less understand where the problem is.

I am not going to give you the full code (yet) because that's not going to do any good. But I will write here the line where I suspect you are running into trouble (straightforward translation to C#):

H(i+N2+1,j+N2+1) = 1./(1+(0.414*(sqrt(((su./Dh)+(sv./Dv)).^(2.*n)))));

I will also suggest an alternative way that you can achieve the same thing, much more straightforwardly, just re-using the stuff I have been talking about so far. Less lines of code, less bugs to worry about.

Now, you might notice that the authors in their paper, state that $H$ should be made symmetric around the origin along the orientation you are trying to filter. So, at the end of the code that is building $H$ you have an extra step (which I suspect you are not doing at the moment) to mirror it. If you try to mirror this function AFTER it has been captured in some $H$ matrix, it is going to be hard work.

For the following to work, I am assuming that in C# you are using some elementary image processing library.

So, write a function that builds you a $y_7$ from above. Nice, clean, simple, one dimensional bandpass filter in the frequency domain. Once you have that, MIRROR IT by something like $y_8 = [y_7(N_{y_7}:-1:1) y_7]$. The indexing convention is $StartValue:Step:StopValue$, so effectively, if your $y_7$ was "Y7VALUES" your $y_8$ is now "SEULAV7Y7VALUES", notice that we skipped the value of 0. Otherwise, the "Y" would appear twice and that would mess the spectrum.

So, your $y_8$ is now one row of your $H$. To turn it into an "image mask", replicate $N$ copies of $y_8$, one below the other. Once you are done, you have your $H$ which if you were to look at as if it was an image, it looks like two parallel lines, two fat bands, one somewhere to the left of the middle of the image and one to the right. This is your $H$ at zero degrees.

To generate the same $H$ at a different orientation. Simply rotate that initial $H$.

DON'T ROTATE A ROTATED VERSION!!!! (so if you did 45 degrees, don't rotate that another 45 to go to 90), because that brings in distortions due to interpolation. Your image processing library should have elementary rotation capabilities. If it doesn't....good luck, you have the Hough transform coming up in the same paper...

I hope this helps.


Implementing OFDM with audio signals I've read an interesting paper called FingerIO
that transmits audio waves in OFDM,
and I wanted to implement and generate sounds using MATLAB.

(I am quite new to signal processing, so please understand my lack of background knowledge!)

In OFDM, there are 64 subcarriers, each with their own frequency band and holding a symbol.
OFDM applies IFFT on the symbols and converts them into time domain with 64 time domain samples.
The 64 time domain samples are sent to the receiver.
The receiver then applies FFT on the received signal and recovers the 64 transmitted symbols.

Below is the snippet of my code to test the data.
I first created a 64X1 matrix that represents subcarriers in the frequency domain, put data in the subcarriers, and applied IFFT.

for i=1:64
   X_Freq = zeros(1,64);
   X_Freq(i) = 1;
   x_Time = ifft(X_Freq);
end

My anticipation was that the frequency of the signal will increase with the value of i since higher subcarriers have higher frequencies.
However, when I ran the code, the frequency increased with the i during the first half (when i=64/2) of the subcarriers, but it began to decrease after the half point.

So I thought, since an FFT returns a symmetric positive and negative frequencies data, I've added the code below before the IFFT operation to create the positive frequency data.

for i=1:64
   X_Freq = zeros(1,64);
   X_Freq(i) = 1;

   X_Freq = [X_Freq X_Freq(end:-1:1)];

   x_Time = ifft(X_Freq);
end

After revision, the frequency now increases with the value of i, but the problem is that the output of the IFFT is now 128 samples.

Does anyone know how I can use 64 samples and use up to 64 subcarriers?
I would appreciate any comments/advises!
Thanks in advance!

You forgot to mix or heterodyn the complex FFT output up to the carrier frequency.  Then the negative index FFT bins will be below the carrier in frequency.  FingerIO appears to use a very high audio frequency carrier.


Does Forward Error Correction (FEC) affect upsampling factor? Is the redundant data introduced during forward error correction (FEC) counted as upsampling?

I'm taking a binary data stream, doing convolutional encoding of rate 1/2 on it, then making it NRZ signal and then upsampling it. In this page, the relation between $EsN0$ and SNR uses a two variables $T_{sym}$ and $T_{samp}$ as $T_{sym}~/~T_{samp}$ which indicates the upsampling factor. So should I include the rate 1/2 of FEC in $T_{sym}~/~T_{samp}$  calculation or not?

Is the redundant data introduced during forward error correction (FEC) counted as upsampling?

Yes.

BUT, under certain conditions.

The whole point of FEC is to be able to detect and correct an error instead of simply detecting it and asking for re-transmission.

To achieve this, we encode the source symbols in a predictable and reversible way, adding some redundant information that, upon reception, helps us infer the symbol that was transmitted even if what was received was incomplete or "damaged" (but always within reason).

This does not come for free. The code rate expresses the amount of redundancy inserted in the data stream. For example a convolutional 1/3 code adds two redundant code bits for every source bit. Consequently, we now have to send 3 bits down the same line, instead of 1. Therefore, our effective data rate is lowered.

One of the simplest ways to "encode" the source data is to repeat them a number of times. Instead of sending 1, send 11111 (a 1/5 code).

On the receiving side, you decode by majority decoding. If most of the symbols received were 1, it is inferred that what was intended to be transmitted was 1 (and conversely for 0).

Upsampling is composed of two steps: The actual upsampling followed by interpolation which is required to "make up" the values of those in-between samples whose value we don't know because it was never captured.

So:

if you were to take a data stream, upsample it and use nearest neighbor interpolation (effectively repeating each sample a number of times) and...
if you were then to send this augmented data stream down a communication channel, and...
if you were to now use a majority decoder on the other side of the channel to infer which one was the level that was actually sent...

then, yes, THAT could count as upsampling.

Otherwise, FEC and upsampling / interpolation are only somehow related as both operations can modify the data rate and content of the signal in a predictable and reversible way. But they do not overlap so much as to count as the same thing.

Hope this helps.


Discrete time non-linear time invariant system dynamics descriptions (state-space or input-output relationship) For the sake of simplicity the following notation $a_k := a[k]$ is assumed for time sequences.

A completely general discrete-time (DT) non-linear(NL) time-invariant (TI) dynamical system can be described with a system of vector difference equations, where $x[k]$ is the system state, $u[k]$ is the system input and $y[k]$ is the system output:
$$\begin{cases} x[k+1] = f(x[k],u[k])\\ y[k] = g(x[k], u[k]) \end{cases} \quad\forall k \in \mathbb{Z} \ge0$$
$f(\cdot)$ and $g(\cdot)$ are general nonlinear functions.

The first equation is very close to a general autonomous non linear differential equation with $\dot x(t)=f(x(t),u(t))$ where $u(t)$ is the source term.

Now, the same dynamical DT NL TI system could be described also through an input output relationship, involving more among their "samples" (current and past):
$$y[k] = h(y[k-1],y[k-2],...,y[k-n],u[k],u[k-1],...,u[k-m]), \quad\forall k \in \mathbb{Z} \ge 0 \quad$$ where again, $h$ is a general non linear function, and $n,m$ are integer positive values.

How can one prove that this input-output relationship form is equivalent to the former state-space representation, at least when $x[0]=0$ (without any care about a possible relationship between $f,g$ and $h$)?

Well, any input-output representation obviously admits a state-sapce form. for your equation in $y[k]$ you can easily construct one as follows. Create a "shift" system (an integrator chain) as
$$
\begin{aligned}
x_1[k+1] &= x_2[k],\\
x_2[k+1] &= x_3[k],\\
&\vdots\\
x_n[k+1] &= y[k]
\end{aligned}
$$
In this way indeed you have $x_n[k] = y[k-1]$, $x_{n-1}[k]=y[k-2]$, ..., $x_1[k] = y[k-n]$. You can do the same by defining another shift system to carry the information in $u[k]$, namely
$$
\begin{aligned}
z_1[k+1] &= z_2[k],\\
z_2[k+1] &= z_3[k],\\
&\vdots\\
z_m[k+1] &= u[k]
\end{aligned}
$$
Thus, with $\xi := col(x_1,...,n_n,z_1,...,z_m)$ the IO representation is equivalent to
$$
\begin{aligned}
\xi[k+1] &= A\xi[k],\\
y[k] &= h(\xi_1[k],\dots,\xi_n[k],u[k],\xi_1[k],\dots,\xi_m[k])
\end{aligned}
$$
The converse is instead not trivial. Pick the LTI system
$$
\begin{aligned}
x[k+1] &= \begin{bmatrix}2 & -1\\0 & 0.1\end{bmatrix}x[k] + \begin{bmatrix}1\\1\end{bmatrix}u[k],\qquad (1)\\
y[k] &= x_2[k]
\end{aligned}
$$
Then the output is the same as
$$
y[k] = 0.1 y[k-1] + u[k],\qquad (2)
$$
however, (2) and (1) are not equivalent.
In general, given a IO representation you can always find a state-space representation, however for general state-space forms a IO representation might be not enough to describe the whole state information. You should find a diffeomorphism among the two representations to assess equivalence


Difference between delta modulation and differential PCM Delta modulation encodes the difference between a sample and the previous sample in only one bit.
Differential PCM encodes the difference between a sample and its predicted value. Is this difference encoded in one bit or more than one bit? If encoding is done in more than one bit, what is the advantage of differential PCM over delta modulation?

Since the encoded values are differences between original values, the values are relatively small values and they are in general distributed as small positive and negative values around zero. This helps to further compress the data by using some algorithm that can use shorter codes to represent smaller values that appear frequently in the data, and longer codes for larger values that appear less often.


Calculating the SNR of Audio Signal (Recommended Libraries) Are there any open source packages or libraries available which can be useful in calculating the SNR(signal to noise ratio) of an audio signal.

The input will be just an audio signal and I have to calculate the SNR of that signal. Can anyone suggest any library in python/C/C++ or any other language which does this or can be useful?

It is really hard to create something to estimate the SNR in all cases. Each case requires its own treatment according to model of the acquisition of the signal.

For instance, if the model is AWGN noise an the audio is of human voice and data is samples in high sample rate (Let's say above 44.1 [KHz]) then you can use a lot of the bins in the DFT of the signal to estimate the Noise STD.
Then according to the energy in the bins of the data you can estimate the SNR.

Regarding a library, you might find general tools in Mozilla's Deep Speech library.


What is the HOMING operator? I was reading "Optical Recognition Systems" book and in chapter 3.3.2.3 where methods for full genetic programming are described, a paper is describing an ADF (automatically defined functions) like this:

  The ADFs have the terminal set Tf ={X, N, S, W, E, NE, NW, SE, SW} and the
  function set Ff = {AND, OR, NOT}. The main program has the terminal set Tc = {I,
  L, NIL} and a set of functions (Fc) containing movement functions, logical operators,
  four ADFs and the HOMING operator.

I am confused about the meaning of HOMING operator. An explanation or some links to point me in the right direction would be great.

The 1993 paper by John R. Koza (a pioneer of genetic programming), Simultaneous Discovery of Detectors and a Way of Using the Detectors via Genetic Programming, can be accessed on CiteSeer. The HOMING operator is described as:

  The one-argument HOMING operator evaluates its argument and, in
  addition, has the side effect of rubber-banding the turtle to its
  position at the start of the evaluation of the HOMING. HOMING is
  equivalent to the brackets in a Lindenmayer system.

From my limited understanding, it is a formal programming operator, within which a certain number of instructions are performed (like motion to North, East, South-West, etc.). It memorizes some initial state (or bookmark), and returns to it when all inner instructions have been executed. In the paper context, it refers to a turtle moving from pixel to pixel to detect a letter in a pixel grid. At the end of its journey, the turtle shall go back to its original location (rubber-banding), or HOME.

Maybe the illustration of  the bracket effect in an L-system (Lindenmayer system) given in section Turtle branching from Drawing simple generative organics with L-systems can help a bit:


ZOH non ideal sampling loss formula Hello i am trying to calculate the amplitude loss of non ideal sampling.
My signals amplitude is 10 at 70 frequency.
When we sample it at Fs=400 we have al 0.511 loss as shown in plot bellow.
How can we get this 0.511 using the ZOH sinc formula bellow. What is f,Tpm,Ts is here so i could see the 0.511 loss? MATLAB code bellow.
Thanks.

f3=70;
Fs=400; % sampling frequency is a bit above 2 times to get all the peaks.
Ts=1/Fs;
Tn=0:Ts:1;
fft_L=length(Tn);
y4_samples=10*sin(2*pi*f3*Tn);
%stem(Tn_new,y4_samples);
ff=fft(y4_samples);
ff1 = abs(ff/fft_L);% normalised FFT
fft2 = ff1(1:floor(fft_L/2)+1); %first half of the vector
fft2(2:end) = 2*fft2(2:end); %we multiply the amplitude by 2
f = Fs*(0:fft_L/2)/fft_L;%frequency vector
plot(f, fft2)

Define this

  Tn=0:Ts:1;

to be

  Tn=0:Ts:1 -Ts;

By defining it the way you have in the program you have not included an integer number of cycles in the samples. This will result in spectral leakage because your sinuoid now does not fit into one frqeuency bin and thus spill out into other bins and this reduces in magnitude.


Partitioning the frequency axis with shifting resolution I am calculating the sum of the signal energy in different frequency bands for music in Matlab. At the moment I have calculated the FFT for a music signal and then split this frequency vector in N partitions. This gives me the sum of the energy for N parts which is all good.

The problem with this is that music most often is concentrated around some frequencies and since my partitions are of uniform length this does not give a good representation of the original signal.

What I would like to do is to be able to give a frequency value f0 and the number of partitions N. This should then give me N partitions but where the partition lengths are shorther around f0 and increasingly longer when farther away from f0. How could one do this?

My idea is that the end result would be a normalized length vector.

Any help would be appreciated.

After some thinking I managed to create an answer that works for my problem. It does not center around the frequency f0 automatically but there is 4 design parameters that you can tinker with. My plan is to run the design parameters through an optimizer. So I can find the parameters that minimize the difference between the sum of the energy in each partition.

N_partition = 20;                                   % number of partitions
N_switch = -2;                                      % switches how how many partitions there are in
                                                    % the lower vs the upper part of the signal
N_lower = N_partition/2 + N_switch;                 % number of partitions in the lower part of the signal
N_upper = N_partition/2 - N_switch;                 % number of partitions in the upper part of the signal
diff_lower_upper = 0.05;                            % design value that decides how big the first partition in the signal is
a_upper = 25;                                       % design value that decides the slope appearance for the lower part of the signal
a_lower = 2;                                        % design value that decides the slope appearance for the upper part of the signal
vec_u = 0:1/N_upper:1;                              % vector for the lower part
vec_l = 0:1/N_lower:1;                              % vector for the upper part

part_upper = ((a_upper.^vec_u-1)/(a_upper-1));      % the change vector for the upper part of the vector, exponential
part_lower = ((a_lower.^vec_l-1)/(a_lower-1));      % the change vector for the lower part of the vector, exponential
part_lower = fliplr(part_lower).*diff_lower_upper;  % the change vector but corrected for the diff_lower_upper design value

change_vec = [part_lower(1:end-1) part_upper(2:end)]; % the change vector, removes the zero value from part_upper and part_lower

nl = change_vec./sum(change_vec);                   % normalize the change vector, this is our normalized length vector

For my case, this creates the following accumulated frequency/partition. I wanted high "resolution" around my f0 that is 1000Hz and then the partitions should grow in size the farther away they come from f0 which they do. Thanks to everyone that commented!


If my sampling rate is not high enough to capture all frequencies, can I accurately capture low frequencies? Let us say that my sampling rate is 1000 Hz. This enables me to accurately capture every frequency in the 0-500 Hz range according to the Sampling Theorem using a DFT. However, there are higher frequencies present in the signal in between 1500-2500 Hz. If I am to convert the 1000 Hz data to its DFT, I understand I will not get the higher frequencies. However, I only care about the lower frequencies. So my question is, will I accurately get the lower frequencies of the signal and their precise amplitudes using the DFT?

If there are higher frequencies present in the signal you will indeed get all these! The problem is the frequency components will be aliased as additional signals also within your first Nyquist zone of 0-500 Hz (assuming a real signal, if complex then the full 0 to 1000 Hz spectrum in the DFT will be unique).

To uniquely observe the signals that are in the 0-500 Hz frequency range, you MUST have an anti-alias filter in the analog domain prior to the sampling of the signal at 1000 Hz. Ideally (but not achievable) you would use a "brick-wall" filter to pass the DC-500Hz band with no distortion or loss, and then completely reject everything above 500 Hz.  Welcome to the real world, where filters have finite impulse responses and we don't have forever to wait for long filter delays-- meaning we must have a transition band in the design of these filters. The amount of effort and cost you want to put into the anti-alias filter dictates how much of the DC-500 Hz band you can reliably observe, and really goes into a lot of consideration in the overall system design (trading design detail elsewhere). A typical number can be 35-40% (35% Meaning DC-350 Hz passband, 650 Hz + stop band.).

Now if your goal is to instead observe the spectrum that is present in the 1500 - 2500 Hz frequency band, you could instead design band pass anti-alias solutions to select them specifically prior to sampling with the same 1000 Hz sampling clock (this is called "under-sampling). In this case you have restricted what will alias into the DC-500 Hz band, and can therefore uniquely identify the actual frequency components from mapping the aliasing. Like above, if it is a real signal and you are restricted to the 1000 Hz sampling rate, then 350 MHz bandwidths would be reasonable for the bandpass filter, meaning you would not be able to observe with one DFT the full 1500-2500 Hz analog spectrum with a real sampling rate of 1000 Hz.

Undersampling is best demonstrated with the graphic below showing two analog signals, an 11 Hz sine wave and a 1 Hz sine wave both sampled with a 10 Hz sampling clock. Once sampled we will only have the red circles representing the signal. This shows that the input in this case could have been either the 11 Hz or the 1 Hz signal, and if no filtering was present to restrict one from being possible, we would have no way to distinguish which signal was present. If both are present then they could destructively or constructively interfere depending on their relative phase.


How to generate IQ samples of a signal? I am using a HW with LMS7002 RF transceiver chip. This chip takes IQ values as an input.
 1. Can anyone point me to a sample IQ value generator? That I can use to test.
 2. What will be the IQ values of a pure sine wave say (1kHz) and how to extract it, say if I try to generate it in software like python or C++?
3. How can I generate a signal that occupies a certain bandwidth say 30 MHz?

The reason one uses a complex signal at baseband is so that one can distinguish between positive and negative frequency...otherwise, when the signal is mixed to a higher frequency, it will have two identical sidebands, taking up twice the spectrum.  Without knowing what you're modulating, I can only offer general advice...take your real baseband signal taking up bandwidth 0..Fmax, and multiply it by exp(2*piiFmax/2*t) (note the 2s cancel out, but the form is more recognizable).  This will heterodyne it down to a complex signal occupying -Fmax/2..Fmax/2, and that can be passed to the RF chip.


DCT: Alternatives to quantization When using the discrete cosine transform are there commonly used alternatives to quantization to decide which dct components to keep/are important? If not, how do people come up with quantization matrices (such as in JPEG), is it purely arbitrary, based on empirical properties of common images, maths, etc.

From what I read there are certain quantisation matrices for different applications, Adobe Photoshop has like 15 or something. The idea is that although there are error based mathematical calculations for these matrices the best way to deduce a good quantisation matrix is by simply using a human eye (in the case of JPEGs)I haven't seen any quantisation matrices that don't follow the general trend of having heavy penalties(High number) on the lower right frequencies, although some images have a large high frequency coefficients the lower frequencies are generally larger.


What is a "jump response" of a filter? I know how to determine if a filter is FIR or IIR, but I don't understand the question. What is a "jump response"?

It's a bad translation of step response, i.e., the filter's response to a unit step $u[n]$ at the input. Note that the step response is the convolution of the unit step with the filter's impulse response:
$$s[n]=\big(u\star h\big)[n]=\sum_{k=0}^nh[k],\qquad n\ge 0\tag{1}$$
where I've assumed causality, i.e., $h[n]=0$ for $n<0$.
From the given step response and Eq. $(1)$ it should be easy to derive the impulse response $h[n]$. You can also immediately draw a conclusion about its length. Then, by considering the symmetry of $h[n]$, I'm sure you know how to derive the group delay of the filter.


linear-expansion filter for image (2D array) I am working in python analyzing 2D arrays and I am looking for a filter or some other routine that would provide a linear expansion of an region boundaries in a 2D array with binary regions (i.e. a black/white image).  For example given this star image:

I would like to apply a filter that expands a region's borders with a  radial, linear decrease from a maximal value (i.e. black) to the surrounding value.   Graphically represented as:

So that the resulting image/array is something like:

Simply resize the binarized image to a larger image but do not use any high order interpolation in resizing and instead use 'nearest neighbor'interpolation technique. In MATLAB:

im2 = imresize(image, 10, 'nearest')

part of 10x expanded image:

and a raster of the image:


How to detect real time sound onsets? I'm a Software Engineer Undergraduate. I'm trying to detect drum strokes(beats) in real-time in python. I used amplitude threshold-based onset detect (beat detect) algorithm but it takes other noises also onsets and it omits some actual beats with low amplitude. Can you have suggestions to optimize my code or another way to solve this?

My guess is that the currently popular and new and likely robust way to solve this detection problem is to feed a sequence of audio fingerprints (such as MFCCs) to an RNN machine learning algorithm that was trained on a large wide range of rhythm tracks mixed with increasing levels of realistic background noise.

Feeding audio stream samples directly to a deeper CNN+RNN ML inference engine might work even better by detecting humanly resolvable beat details that audio fingerprinting misses, but that would be computationally far less efficient and would take a lot longer to train.


What does the term "railed" mean in signal processing? I'm having trouble finding a definition of "railed" that relates to signal processing.

Am I correct in my guess that this term is in fact from this field?

My signal data comes from an EEG device. The lightly documented open source software I'm using doesn't define it, but it shows that term when there is no signal data being displayed.

Is that all it means (no data)? Or does it mean something like the signal being read is too great to be displayed or correctly measured?

A railed signal, or a railing signal, seems to indicate a flatline. On BIOPAC, Railing signal (flatline) says:

  When the amplified signal for any given channel exceeds the range -10
  to +10 volts, the signal will rail. You will see a straight line at
  -10 or + 10 volts (more likely the reading will be close 9.99 volts). The MP system is designed to work only in the range -10 to +10 volts.
  The signal could rail for several reasons (which are not exclusive)...

From Amplifiers: What do rail-to-rail and single supply mean?

  With respect to analog signals, a “rail” is a boundary that a signal
  has to work within.


Speech synthesis requiring very little CPU performance? Back in the days of 1 MHz 8-bit CPU personal computers (Apple II, Atari 800, et.al.), there were software programs that could do comprehendible arbitrary text-to-speech synthesis on those PCs.  What published speech synthesis algorithms might be suitable for implementation on micro-controllers with similarly limited performance and memory?  (If any...)

The methods I played with on 8-bit home computers back in the 80's involved having a small set of phonetic samples and then having the text you want to speak defined in terms of those samples.

Since I didn't know about IPA at the time, I just used two sounds per letter and allocated the short sound to lower case letters and long sounds to the uppercase. This did leave quite a few gaps though, o for /ɒ/ (hot) and O for /oʊ/ (open) were easy, but /uː/ (the oo sound in moon) couldn't be represented in my system, which was a problem for me, given my name (/buːð/ rhymes with /smuːð/, try them in itinerariums phoneme synthesis web page).

These days, you could do much more clever stuff. For a start you could use lookup techniques to translate known words in orthographic text into your known phonetic symbols, and for unknown words, implementing Phonics as an algorithm should do a reasonable job.

After getting your phonetic stream, you could then probably use DMA techniques to stream individual phoneme samples from memory to your DAC, with less overhead to your CPU.


What steps come after finding the gradients of the SIFT key points? I am using a SIFT algorithm to extract features from an image.  I understand that the SIFT descriptor first finds the extrema and then finds the gradient and direction of each of these interest points.

What happens after this part?

First of all, there are two distinct parts to SIFT.  The first part is interest point detection algorithm (aka key-point detection), which finds local extrema of the multi-scale difference-of-gaussians function.

The second part is computing the feature descriptor, which is a vector describing the image patch around each key point. SIFT computes this descriptor by taking a 16x16 block of pixels centered at a key point, dividing it into 4x4 cells, and computing an 8-bin histogram of gradient orientations within each cell.  This results in a 4*4*8=128-element vector, which is the descriptor.

Finally, the SIFT algorithm normalizes the feature descriptor for rotation.  It finds the dominant gradient direction in the patch around the key point and shfits the gradient histograms such that the first orientation bin corresponds to that direction. This makes the descriptors invariant to in-plane rotation.


Phase spectrum of a Fourier transform I am trying to compute the phase spectrum of the signal
$$
s(t)=\frac{A}{\pi}\left[H\left(t+\frac{\tau}{2}\right)-H\left(t-\frac{\tau}{2}\right)\right]
$$

The Fourier transform is
$$
S(j\omega)=\frac{A\tau}{\pi}\cdot\frac{\sin\left(\omega\displaystyle\frac{\tau}{2}\right)}{\omega\displaystyle\frac{\tau}{2}}$$

My answer for this would be $\phi_s(\omega)=0$, since the imaginary part of $S(j\omega)$ is zero, but according to my textbook, the answer is

$$\phi_s(\omega)=\begin{cases}0&S(j\omega)>0\\\pm\pi&S(j\omega)<0\end{cases}$$

How can I get to this answer and interpret it correctly?

A real-valued negative number has a phase of $\pm\pi$, because $e^{\pm j\pi}=-1$. So for all frequencies for which the real-valued spectrum is negative, the corresponding phase angle equals $\pm \pi$:

$$S(j\omega)=\begin{cases}|S(j\omega)|,&&S(j\omega)\ge 0\\
|S(j\omega)|e^{\pm j\pi},&&S(j\omega)<0\end{cases}$$


How to calculate this energy for a given signal? So the question Im asking can be seen in the picture below.The d at the end of the question stands for the minimum euclidian distance.

Mostly I just want to know how to calculate the energy for each signal. I know how to calculate the energy efficiency between two constellations (10 * log(E1/E2)).

Thank you

If your question really is how to calculate the energy of a given signal $x(t)$ then the answer is

$$E_x=\int_{-\infty}^{\infty}|x(t)|^2dt\tag{1}$$

but I wonder if that is actually what you're asking, because it would have been very easy to find that formula.


2 points and instantaneous phase sufficient to estimate amplitude and frequency of sinusoid? If you have 2 points (x1y1,x2y2) in a Cartesian grid, there are an infinite amount of sinusoids that pass through them. But if you are also given the instantaneous phase (Th0) in radians or degrees at that point, is that sufficient to estimate the amplitude and frequency of the sinusoid that passes through these 2 points? And is there a way to find a formula f(x1y1,x2y2,Th0) where you can simply plug in these 3 variables which returns the amplitude and phase of the sinusoid, so you don't have to look it up in a table? Thanks.

An additional requirement is that the two points not be aliased in terms of their equivalent sampling rate relative to the frequency of the sinusoid.  The distance between the two points can't be a equal or greater than half the period of the sinusoid, or the solution could be ambiguous.


Jigsaw puzzle: isolating the pieces (separating "stuck" blobs) I'm new to image processing, and I'm trying to get my feet wet. I have taken a picture of some jigsaw puzzle pieces and I want to isolate the pieces.

I am using the Python SimpleCV library to do this, and so far I've managed to get pretty decent results using findBlobs(), hueDistance(), and drawMinRect(). Here's what I've got so far

Pretty good, honestly. The only place it's wrong is with those two pieces touching slightly to the right of center.

I tried dilating the image, but that seems to exacerbate the problem. How can I "shrink" these blobs to eliminate the overlap?

I agree with Abid Rahman K, however eroding only would result in shape shrinkage.  To only get rid of the gap, Abid Rahman suggests a watershed based approach, but for a faster and easier method you might try:

Erode the regions
Apply connected component labeling.
Dilate the regions back.

This way you will have almost the same shape (with some artifacts arising due to morphology) but you will surely end up with single connected components.

Watersheds and distance transform are also powerful and might be good choices if you have a lot of overlap.


Design FIR with linear phase response What is the best way to design a FIR filter with a given linear phase response, such that each frequency inside the passband is phase shifted according to:
$$H(f)=e^{i2\pi f k}$$
where $k$ is  factor that I want to specify. I also want to specify the number of taps $N$. Note that the effective shift is bigger than the tap length $N$.
I don't really care about the bandwidth of the passband, but it has to be at least 1 MHz at a sampling of 100 MHz. The amplitude should be relatively flat inside the passband.
In python there is a lot of functions which allow me to specify a passband and how flat it is, but I have not found any way to control the phase according to my needs.

What you want is basically a delay
$$x(t-\tau) \Rightarrow e^{-j\omega \tau}$$
If you can live with $k$ being quantized you can simply use a single tap delay FIR filter. If you need more granularity you need to implement a fractional dealy. I would cascade a shot fractional FIR filter with a long single tab bulk delay filter.
If $k$ is positive the filter becomes a non-causal. The best way to deal with this depends a bit on your application. If you do off-line processing, you just need to delay all other signals paths by the same amount.


Echo effect on sinusoidal input: Intuition I am currently reading about adding echo effects on signals. I was wondering if my understanding is correct.

In the attached figure, the blue line represents the original signal and the red line is the delayed original signal. Is it correct to add these two signals to represent the echo effect, resulting in the green line, assuming no attenuation?

I understand that an echo effect also attenuates the signal. However I just assumed that there is no attenuation for simplicity.

Is it correct to add these two signals to represent the echo effect, resulting in the green line, assuming no attenuation?

Yes

  I understand that an echo effect also attenuates the signal. However I just assumed that there is no attenuation for simplicity.

To continue along your example, if you delay the signal half the period of your sinusoid, the "memorised" positive front from the delay will be added with the negative front of the second half of its period and the result will be zero.

This exact mechanism is responsible for the frequency response of the comb filter.

You can create more complex environments by using an impulse response and convolution.

Hope this helps.


Pitch perception - Why does our ear not directly detect the missing fundamental frequency? My motivation for this question is related to biological phenomenon of our inner ear, so I initially posted in the biology community, but I suspect the answer has a DSP explanation, so I figured this might be a better place to get an answer.

I'm learning about pitch perception in an Intro to Music class right now, and learned about the case of the missing fundamental.

In the main image in that wikipedia page, it seems like in the bottom graph, with the fundamental frequency and its second harmonic removed, the wave is still very periodic at 100 Hz.  Since the Organ of Corti has some area which should be excited when there's a sound wave of frequency 100 Hz, why does this soundwave, which seems periodic at 100 Hz, not directly excite that region of the organ?

(like, won't there be high air pressure every 0.01 seconds, which will sway that organ in some direction, and thus excite that region which is "resonant" for sounds at that frequency?)

Maybe another way to phrase the question: why does the top signal have a signal with frequency of 100, but the bottom signal does not have a part with frequency 100?  It seems to me that the bottom signal also has a part with frequency of 100 Hz.

Here is roughly how it works.

A typically musical note is made out of harmonics and the fundamental. The fundamental is also the spacing between the harmonics.
For 100 Hz fundamental you would get 100Hz, 200Hz, 300Hz, 400Hz, 500 Hz, etc.
For 200 Hz you would get 200Hz, 400Hz, 600 Hz, which is simlar to 100Hz but not the same. Note that specifically 300Hz and 500Hz are missing. The spacing here is 200 Hz
In practice you  rarely get all the harmonics, some are louder than others, some are missing all together, however, you still can easily hear the pitch
That's because you determine pitch by the spacing of the harmonics, not by the lowest.
For example, if you generate a sound by adding 500, 600, 700 &  800Hz, you will hear the pitch corresponding to 100Hz, although no frequencies below 500 Hz are in the signal at all. 100 Hz is the periodicity, not necessarily the frequency content.
The organ of corti does actually a mechanical Fourier transform. The excitation through the oval window creates a bending wave on the basilar membrane. Due to mechanical resonance of the membrane and some active feedback mechanism controlled through efferent nerves and hair cells, there is a strong localized peak in the bending wave, the position of which corresponds with the frequency
In our example above you would see peaks at 500, 600, 700 & 800 Hz positions, but NOT at 100 Hz. There is no 100 Hz in the signal, so there is no excitation at the 100 Hz location.
The actual pitch detection is done by identifying correlated spectral parts parts of the signal and determining if there is a constant spacing. The spacing is the pitch. I don't know of the top of my head where in the periphery that happens but I'm guessing it's most like down through some correlation mechanism


How do you interpolate between points in an image (2D), e.g. using splines? I can understand just fine how to use 1-dimensional interpolation on data points where one coordinate is a function of the other: y = f(x). However, when we have an image, we generally have points in the image that are not given uniquely by such a function, e.g. two points (x1,y1) and (x1,y2), so how do we describe this in order to interpolate between points in the image?

There are various methods on 2d interpolation (this one, and this one). But most of them considered at least 4 points rather than 2. The simplest 2d interpolation is 3 1d interpolation, in which you interpolate the points between (x1-d1, y1-d2) and (x1+d1, y1-d2) as (x2,y2), then you interpolate the points between (x1-d3, y1+d4) and (x1+d3, y1+d4) as (x3,y3). Finally (x1,y1) is obtained by 1d interpolation of (x2,y2) and (x3,y3).


windowed sinc filter in matlab Hello i'm designing a low pass filter windowed sinc in matlab

Fs = 2048;
Fc = 250/Fs;  % 250hz
M=2048 % (to get a 4hz transition bandwidth) sig contains 200hz + 300hz sinusoids, i'm aiming at leaving only the 200hz

filt =sinc(2*Fc*((1:length(t))- M/2)).*(0.42-0.5*cos(2*pi*(1:length(t))/M)+0.08*cos(4*pi*(1:length(t))/M))';
filt(M/2)=2*pi*Fc;
fsig = sig .* filt;

i understand that the winsinc filter is applied by convolution in the freq domain, so it means a multiplication in the time domain, which i did, but i get an erroneous result, i get a 0 flat frequency everywhere in my frequency plot.
what's wrong?

doing a freqz(filt) i get this, it seems correct reponse

Your understanding is mistaken. LTI filters are applied via convolution in time domain and multiplication in frequency domain.

The confusion possibly happens about designing the filter vs applying it.

The windowed sinc filter impulse response is given by

$$ h[n] = w[n] h_i[n] ,$$

where $h_i[n]$ is the ideal (infinite length) sinc filter which corresponds to a brickwall lowpass type in your case. To approximate and taper that filter, you apply a window $w[n]$ of length $L$ to it by multiplication in time domain and its effect in the frequency domain is a convolution by the well known Fourier modulation theorem:

$$H(e^{j\omega}) =  W(e^{j\omega}) \star H_i(e^{j\omega}) $$

However, when the filter is applied on an input, the output equation in the time domain is

$$y[n] = h[n] \star x[n] ,$$

and the corresponding frequency domain equation is

$$Y(e^{j\omega}) = H(e^{j\omega}) X(e^{j\omega}).$$

So, rewrite your code base don this to get an expected result.


Extracting common signal without knowledge about noise Given two noisy time series thought to contain a common signal,
$$ x_1(t) = s(t) + n_1(t), \quad x_2(t) = s(t) + n_2(t), $$
what is the best way to determine $s(t)$ without assuming any distribution for the noise terms $n_i(t)$?

In particular, the $n_i(t)$ are neither Gaussian nor stationary. The only assumption I am willing to make is that $n_1(t)$ and $n_2(t)$ have zero mean and are uncorrelated with each other, i.e. $E(n_1(t) n_2(t)) = 0$ (though even this assumption can be broken in my problem, but only slightly).

Is it true that the best estimate is simply
$$ \frac{x_1(t) + x_2(t)}{2}? $$

I am particularly interested in finding not only an estimate for $s(t)$, but also a measure of the uncertainty of this estimate.

Taking the average is a reasonable heuristic. The term “best” needs a lot more clarification because it needs a criteria.  Assuming that the signal is exactly the same on both channels is actually a strong assumption, so if you consider small perturbations to your model, say a time delay between the two channels, which is a lot more reasonable than perfect alignment, the character of the problem has changed.

There are techniques like ICA but they work well when they do and don’t when they don’t.

If your question is of a theoretical nature, there are answers to questions like is there a common signal in both channels using non parametric statistics.

If this is a practical question, then you should look at your data very closely and see whatever else you can reasonably assume.


Zero Padding in Frequency Domain (Spectrum): How to Split the Middle Bin (Nyquist Frequency)? When a complex digital signal is converted to its spectrum through an FFT, the result will contain a series op positive and negative frequencies. That is, when the spectrum has N bins. Bin 0 will be the DC, bin 1->N/2 will be the positive frequencies and bin N-1 -> N-N/2 will be the negative frequencies. In this case it is clear that bin N/2 represents the radial frequencies $\pi$ and $-\pi$ at the same time.

Now, If we zeropad the middle of the spectrum (e.g: when upsampling) such that its new size will be M, we obtain a new one as follows. $X$ is the small spectrum. $Y$ is the zero padded larger spectrum.

$X_{[0:N/2[}$ is mapped to $Y_{[0:N/2[}$
$X_{]N-N/2:N[}$ is mapped to $Y_{]M-N/2:M[}$

In the above, I skipped the middle bin ($X_{N/2}$) because it is not yet clear whether it should assigned to the positive or negative section of the new signal.

To better understand what would be the best result, I performed a small test by creating a random spectrum of length M, then converting it to the timedomain, picking out every M/N-th sample and then converting it back to a spectrum, but this time of length N. The result showed that the following is true for any random spectrum:

$X_{N/2}=Y_{N/2}+Y_{M-N/2}$

Further testing revealed that indeed, we can split the value of the middlebin however we want, as long as their sums match the original value. Yet, the interpolations between the samples (when converting $Y$ back to the timedomain) oscillates differently for each particular choice.

The following three images show the effect of the middle bin on the phase of the upsampled signal. The first one places the middle-bin at N/2 in the new spectrum. The second image splits the content over bin N/2 and M-N/2 and the last one maps the middle bin solely to M-N/2

As you see the interpolation goes each time exactly through each of the blue values. The only difference is how the interpolation oscilates.

Is there any rule-of-thumb (and why), or any argument/standard practice on how to split the value of the middle bin when zero padding ?

One potential issue here is that you are technically in violation of the sampling theorem. Having non-trivial energy at Nyquist is typically an indication that the signal wasn't sufficiently low pass filtered before being sampled. So chances are you already have aliasing. If the Nyquist energy is low enough, it doesn't matter how you split it.

Zero padding past Nyquist to upsample introduces a frequency domain "brick wall" filter between Nyquist and Nyquist + 1. You can certainly go through with the mathematical procedure but I would not expect this to closely match the behavior of the original signal.


Simulate a simple multipath channel for OFDM Modulation test. MATLAB Hi I need to simulate a simple multipath channel in matlab through which I can send my OFDM transmission signal, I want to test my OFDM receiver design.
So far I have:

%% Passing through Channel
nIFFT = 4096;

% Create channel
channel = randn(1,nIFFT) + sqrt(-1)*randn(1,nIFFT);

% cext_data is the signal to be transmitted
after_channel = filter(channel, 1, cext_data);
wavwrite(after_channel,'ofdm_tx8khz.wav');

I get a warning when I run the transmitter:

  Warning: Data clipped during write to file:ofdm_tx8khz.wav

and when I run the receiver I get:

  Index exceeds matrix dimensions.

  Error in ofdm_rx (line 26)

  chunk2 = signal((s+nIFFT):(s+nIFFT+nCyclic-1));

I was really hoping someone could help me simulate a multipath channel for my OFDM program.

Wav files expect values from -1 to +1.  You are probably exceeding those limits when you filter the signal through the channel model.  After the channel filter add the following:

after_channel = after_channel / max(abs(after_channel));

I suspect, though, that you will still have a problem.  I believe that wavwrite ignores the imaginary part of complex data, so you will be losing that part of the signal.


How the select the representative frame when detecting vehicles I am detecting vehicles from the video (taken with a camera), and the vehicle detection part works fine.

My problem is that I should always extract only one representative frame per vehicle, and not the whole sequence.

As an example, if I have a few seconds of video (and possibly multiple vehicles), my approach currently extracts multiple frames (for every vehicle). Instead, as an output I would like one representative frame in which the vehicle is whole, i.e. well visible.

I have already implemented an entropy technique on it which makes it better but am still getting too much frames of the same vehicle. Could you suggest a technique which I can use to extract only one frame which contains the whole vehicle and ignore all other frames which contain that (same) vehicle.

I use the following entropy technique in my problem

double entropy = 0.0;
for (int i=0; i<histNorm.rows; i++)
{
    float binEntry = histNorm.at<float>(i,0);
    if (binEntry != 0.0)
        entropy -= binEntry * log2(binEntry);
}

I can not give you a full approach, but I might have an interesting suggestion. As soon as you detect a vehicle, start tracking it.

For every detection, you could "remember" the first and last frame of the video where you tracked that vehicle. You associate every vehicle detection with a series of consecutive frames. The beginning frame would be the one when the detection was made, and the ending would be the first frame when the tracking algorithm fails.

Then, you could select only a single frame for every distinct vehicle you tracked during your video sequence, with more or less simple approaches:

select the middle frame of the tracking sequence, as there is a high chance that in the middle of the tracking, the vehicle was well visible
do matching between several points in the sequence, or just base your decision on the number of features you can still track at a certain frame:  selecting a frame with a large number of trackable/matchable features is probably a good choice
remember the "quality" of your detection for every frame in the sequence, and select the frame where the detection quality was the best
depending on the preferences, be mindful to either select one representative frame if the sequences for two vehicles overlap, or purposefully select different frames (possibly from non-overlapping parts of the sequences).


Should I use ideal or non-ideal filters for offline filtering? I've got an offline signal that I want to high-pass filter. Should I use a Butterworth filter, or could I use the fact that the whole signal is known and use an ideal (step) filter?

Using a high-degree Butterworth filter is, I think, practically the same thing as using an ideal filter, and a higher degree means a better filter if there's no implementation and processing cost. Could the sharp change in impulse response cause a problem, though?

The signal is a bioelectric signal that is almost periodic, and I've got 100 periods or so.

For most practical applications, filter responses that approach the behavior of the ideal "brick-wall" filter are overkill. I know it's tempting to try to design a really sharp filter when you've got all the time in the world (i.e. for offline applications), but if you really look at the characteristics of your problem, you can most likely get away with something much more reasonable.

Another good reason not to use a brick-wall filter: its impulse response is infinite in length, and it has the form of a $sinc$ function. When you apply such a filter to your signal, you may notice long-duration ringing in the signal at the output of the filter. This comes from the fact that the $sinc$ function in the filter's impulse response is infinitely long and doesn't decay very fast; the resulting effect is most likely not desirable. In general, filters with long frequency responses are not very suitable for analysis of short signals, as the filter output is dominated by the transient behavior of the filter.

The aforementioned effects aren't limited strictly to ideal filter responses; if you design a filter with a really sharp cutoff, it's still possible that you will get time-domain artifacts that you don't want. This intuitively makes sense, because you can view the frequency response of the "not-quite-ideal" filter as the response of the ideal filter convolved with a narrow main lobe that smears the response out a bit. Frequency-domain convolution is equivalent to windowing in the time domain; if you look at the impulse response of a very sharp filter that you've designed, it's likely to look much like a $sinc$ function with a window applied to taper the impulse response off at a faster rate than the $sinc$ does on its own.

I give the same advice as with most filter-design problems: you should really try to tackle the problem as quantitatively as possible. Using intuition on what might seem like a good approach can often steer you down the wrong path. Instead, think about the following:

Where is my signal of interest in the spectrum?
What unwanted signals are in the spectrum? Where are they?
How large are the unwanted signals relative to the signal of interest? In order to accomplish my ultimate goal, how much must it be suppressed?
How much distortion can I tolerate on my signal of interest (both in amplitude and phase)?
What computational limitations do I have?

The first four questions will give you filter performance specifications that you can use with any number of design methods to arrive at filters that will achieve your goals. The last question is also important, and can be used to choose between different filter topologies (i.e. FIR versus IIR) to find one that is implementable for your application.


Formula for calculating IIR filter coefficients for a high order low pass? I'm trying to make a filter for use in real-time audio processing and I'm trying to figure out how to produce coefficients for a low pass with a steep attenuation curve. I've found a few examples of b0, b1, b2, a1, a2 but I'd like to have the option of a high order filter, which to my knowledge means more coefficients. How are the extras calculated? I'd be extremely grateful for any help on this. Thanks!

First, the equations you solve for the coefficients depend on what type of filter you're after. They also depend on how the equations are derived. For instance, a common case would be a Butterworth lowpass, converted from an analog prototype in the s-domain via the bilinear z transform. Converting a higher order lowpass with the bilinear transformation would yield a solution with more coefficients.

However, it's common practice to form higher order filters from second order sections in series, with a first order section if the filter order is odd. The responses of the individual filters multiply to give the overall higher order response. The reason for reducing a high-order filter to first and second order filter is that high-order filters are more sensitive to the limited precision of typical floating point fixed point math. The cascaded filters are all set to the same corner frequency, and the Q of each filter is set to a different values that, together, produce the intended overall Q.

Here's an article about how to produce higher-order Butterworth filters for first and second order sections, including a calculator that produces the Q values for filter orders up to 20:

Cascading filters

Here's a second-order filter calculator that you can use to generate coefficients for the individual filters:

Biquad calculator v2

Here is more detail on the bilinear z transform:

The bilinear z transform


Can a 1-dimensional signal be in a domain other than the time or frequency domain? In signal processing I have only heard of signals being in the "Time" domain or in the "Frequency" domain, but are there situations where a signal is not in either of these?

If so, what are these other domains called? How could you translate a time or frequency domain signal into one of these other domains? Are there any practical use cases for doing so?

This question is referring specifically to 1-dimensional signals (for example a signal read from an audio file, or a simple sine wave).

Really, there's nothing special about time and frequency domain.

The math doesn't care whether you're transforming amplitude over time, gravel over mountain height, or smell intensity over fridge temperature. To the math, you map elements from a function space to elements from a function space when you do the Fourier transform. Full stop.

Any scalar physical entity can be used as basis for a 1D signal. Spatial signals are the most obvious: what's the height profile of something over distance, or the temperature over distance, or color over distance, whatever over distance. In some contexts, the resulting entities are called "spatial frequencies", but sometimes something else. Names are just names.

Or what is the temperature of something over pressure? That's another 1D signal. And it makes a lot of sense to Fourier transform it, if you're dealing with compressible media in a container, things will oscillate, and not necessarily over time. Do I know whether that has a name? No. Did some PhD give it a name in his dissertation about pressure in containers? Pretty likely, if you ask me. Does that name matter to anything but the notation in that dissertation? Maybe. Maybe not.

Very commonly observed: impulse / state space its Fourier "dual" location space in solid state physics.

We can go on with this all day, because, as said anything can be a 1D signal, and some have Fourier transforms with special names, other don't. It's by pure coincidence that you've heard of the time/frequency pair first. If you were a math student, you might have heard about probability densities and characteristic functions first – yet another 1D signal/FT pair.

So, really, to cite someone I think Tukey,

  One can Fourier transform anything – often meaningfully.

Your question is just asking for an arbitrary list of names. Names don't actually matter – they just make it easier to communicate about something. Whether I call something a "time signal" or a "frequency signal", however, makes little difference. For all we care, I could be calling it "original space" and "image space" (under the Fourier transform), because that's what "time" and "frequency" are. It's just convenient that these two have a name.


How Are Unmeasured Properties (Velocity and Covariance of Velocity) Handled with a Kalman Filter? I'm trying to understand how I can update a Kalman filter with a state variable for position and velocity when I only measure position. I have a covariance matrix of the position measurements. But what does the overall covariance matrix of the measurements look like where velocity is unmeasured?
It seems that zeros would send the gain with respect to velocity to favor the measurement, which would drive the velocity estimates toward whatever value I input for the velocity measurements.
I have a similar problem where I am trying to fuse alternating measurements where I either have X and Y but no Z or X and Z but no Y. I would expect a single filter to fuse the two sets of measurements if I do this right, but I don't have a full grasp on how to handle to covariance terms involving the unmeasured property.

This is exactly where the Dynamic Model comes into play.
The whole idea of the Kalman Filter is that you have a model which connects between variables which are measured to those which are not measured (Estimation) or measured differently (Fusion).
Since the velocity is the derivative of the location over time you have a model which connects them both. Once you set the model (Piece Wise Constant Velocity / Changing Velocity and the level of noise) the Kalman Filter does the fusion of the data and the model to estimate the parameters.


How to validate that an audio algorithm result is independent of microphone device? I have a machine learning algorithm that takes speech sample audio recordings collected from mechanical Turk. During processing it was shown that some audio from certain OS/microphone devices have been preprocessed by noise suppression algorithm by default which affects the results I am getting.
As part of my writeup I need to ensure that OS and Microphone preprocessing is accounted for. Ideally I would like my algorithm to be independent of the recording source, however I am not even sure how to go about detecting the kinds of preprocessing that is being applied and give a lower and upper error bounds for how it affect my algorithm.
My question is how does other people deal with preprocessing on audio files that are outside your control? Are there some common preprocessing that I should know about (echo cancellation, noise suppression, etc.)?
Something I can do for example is to add noise suppression myself on some reference files and add +/- for the max and minimum deviation from the results from the reference file.  However it's a bit unsatisfactory as I do not know how different my preprocessing is to other devices.

There isn't a universal pre-processing chain, but here are some common ones that I'd look out for:
1) Noise Gating
If a signal is below a certain threshold, the mic is muted. To identify: You'll see a VERY quiet noise floor, followed by a sudden spike in level.
2) Automatic Gain Control (AGC)
Mic volume modulates based on the signal level. This will be difficult to identify unless you have multiple recordings from the same system, which will show modulation in the noise floor based on the signal level.
3) Compression
This one should be self-evident. You'll see a very steady signal level that makes the waveform look like a "block", i.e. not much changes in level throughout the entire recording.
These are all relatively simple real-time processing examples. There may also be more sophisticated post-processing applied in some cases, such as noise suppression.


What is temporal neighborhood? I read in a paper http://tcts.fpms.ac.be/publications/regpapers/2011/nrmmbgtd_icvs_LNCS.pdf about spatio-temporal filters, which is applied to frames of the videos. I could not understand these terms in the paper if anyone can help.

3D Low Pass Filtering

Temporal Neighborhood

Spatio Temporal Filters

The temporal neighborhood refers to data nearby in time. That is; data located in the frames just before and just after the current frame.

The spatial neighborhood on the other side refers to the data nearby in space. In this case this is the pixels nearby in the same frame.

Spatio-temporal filters (which is 3D filters) takes both nearby pixels in the current frame and nearby pixels in nearby frames into account. In this way you can take advantage of information from more than one frame, which is often a good idea since most video sequences does not change much from frame to frame.


High resolution spectrograms I have a recorded audio signal that I would like to convert into spectrogram with the highest possible yield of information.
The audio clips I am working with is about 30 ms long and contains frequencies between 20 Hz - 20 kHz. The sampling frequency used is 44.1 kHz.

What I have done:

In Python I have used the library librosa to create amplitude spectrograms. Here I have used the length of the signal as number of points for the FFT, hop length (number audio of frames between STFT columns) of 1 and window length (Each frame of audio is windowed by window()) of 64.

For reference, the code used to generate and display the spectrogram is as follows:

plt.figure()
librosa.display.specshow(librosa.power_to_db(
                         np.abs(librosa.stft(signal, len(signal), 1,64)) ** 2,
                         ref=1.0), y_axis='linear', sr=fs, cmap='coolwarm')
plt.axis("off")
plt.savefig(SAVEDIR_spec_cw+"spec_"+str(i)+".png")
plt.close()

This seemed (visually) to give the best image result as seen below:

Is this as good as it gets?

I do have access to MATLAB as well, if someone can suggest a better method using that instead.

Good for what purpose?  There are always trade-offs.

You can play with the color palette or function, or do local non-linear color selection.

You can try compositing/stacking multiple resolutions of spectrograms (both time and frequency) for higher contrast and perhaps more interesting appearances, but that process can add its own artifacts.


Degrees rotation based on frequency offset? What is the relationship to calculate degrees of rotation in the constellation given a baud rate and frequency offset?

frequency is the derivative of phase over time; phase hence frequency's integral over time. And that's all the knowledge you need, since baudrate gives you the time between consecutive symbols.


FMCW radar: Linear phase change determination with sinusoidal decomposition I am struggling to understand how linear phase shift can be determined with a transform whose main purpose is to determine frequencies rather than rate of change of a function.

In FMCW radar, received signal corresponding to a chirp is passed through a mixer and the resulting beat frequency is sampled by means of ADC. The samples are stored in a row of a matrix. Similarly, samples from  other chirps are stacked up in corresponding rows. Different objects at different locations produce different frequencies (thanks to Doppler). So a FFT of each row yields dominant frequency bins. The frequency bins directly relate to object range. Absolutely fine so far.

Assuming that these objects move (and slowly), they must continue to appear in the same frequency bin (column)for all of the chirps. But since they are moving, the phase of the frequency component must change. For simplicity, we assume that the magnitude of all frequency components in a column is constant.

Thanks to Shift Theorem, a signal delayed/advanced in time results in a phase shift of the frequency components of the time domain signal.

Question-1: Am I right in assuming that a constant velocity results in a LINEAR phase shift of the frequency component across chirps?

Question-2: If the phase shift is linear, how can FFT along a column determine the same? FT/DFT essentially performs decompositions in terms of basis sinusoids. But in our case, we are interested in calculating the slope of the line (or rate of change or derivative) representing the phase-shift change across chirps. Derivative is not same as frequency. I think we need a Non-sinusoidal decomposition hich can yield various slopes. But second-stage FFT along the column seems to work.

What have I misunderstood?

Your understanding is mostly correct but no need of non-sinusoidal decomposition and all.

  Am I right in assuming that a constant velocity results in a LINEAR phase shift of the frequency component across chirps?

Yes, the change in phase due to time change $\Delta \tau$ is $\Delta\phi = 2\pi F (2\Delta \tau)$, $2\Delta\tau$ because signal has to travel back so double the distance covered by object. $F_c$ is the chirp frequency. $\Delta \phi = 4\pi F\frac{\Delta d}{c}$, where $\Delta d$ is the distance increment due to velocity of object and $c$ is velocity of light. Since wavelength $\lambda=F/c$, and $\Delta d=vT_c$ (where $T_c$ is chirp duration), phase change between 2 chirps is   $\Delta \phi = 4\pi \frac{vT_c}{\lambda}$.  So you can see that relative phase of peak bin of FFT output is linearly changing with respect to velocity $v$.

  If the phase shift is linear, how can FFT along a column determine the same?

Once you have identified the column corresponding to peak bin, take the FFT of phases of these peak bins column wise. If phase changes linearly with time, the FFT of phases would show up as a peak bin when you take FFT along the column.

From the previous formula, we have $\Delta \phi$ as the change in phase between two consecutive chirps. So $\Delta \phi$ can be denoted as angular velocity $\omega$. Hence $v = \frac{\lambda \omega}{4\pi T_c}$. So if you take FFT of phases, if you see a peak at $\omega=\omega_0$, the velocity is $\frac{\lambda \omega_0}{4\pi T_c}$. If you see 2 distinct peaks, implies there are 2 equidistant objects approaching or retreating at different velocities at that point of time.


Calculating an image's fourier spectrum by hand? Suppose I have a $4x4$ image with the following values as its grey-level intensity for each pixel like this:

I want to get its Fourier spectrum. Usually, I would just punch into Matlab and run a fft for everything to be calculated automatically. Now, if I don't have Matlab, and I have to do it by hand and with just an ordinary scientific calculator, can I still do it?

I know I could do the brute force way on the Fourier transform's equation like this:

In this case, $N=M=4$.

So, $f(x,y)$ returns the intensity value of the above image at position $x$ and $y$. Eg, $f(1,1)=1$, $f(2, 0)=3$

But this will be insane to work out the summations so many times on paper. Most calculators can't do summations with complex numbers too.

I have tried doing it with Fourier transform's separability property, but doesn't seem to work as I still have to go through the tedious summations.

What other ways can I use to calculate Fourier transform on a simple image without using Matlab, particularly for the values in the first row of the example image above?

Your particular case can be solved without a calculator. Fourier transforms are linear operations. You can decompose your input as:

$$3 \left( \begin{array}{cccc} 0 & 0 & 1 & 1 \\
0 & 0 & 1 & 1 \\
0 & 0 & 1 & 1 \\
0 & 0 & 1 & 1 \end{array} \right) + \left( \begin{array}{cccc} 0 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \end{array} \right) $$

The first matrix is constant on the Y axis, so its 2D FFT is only the 1D-FFT of one row, times 4, and 0 everywhere else. The 1D FFT of $(0, 0, 1, 1)$ is $(2, -1+j, 0, -1-j)$, so this gives you, for the first term:

$$3 \times 4 \times \left( \begin{array}{cccc} 2 & -1+j & 0 & -1-j \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \end{array} \right)$$

As for the second term, there's only one value so there's only one component in the sum - it's just one of the basis functions, the checkerboard pattern:

$$\left( \begin{array}{cccc} 1 & -j & -1 & j \\
-j & -1 & j & 1 \\
-1 & j & 1 & -j \\
j & 1 & -j & -1 \end{array} \right)$$

Summing these two, you get:

$$\left( \begin{array}{cccc} 25 & -12+11j & -1 & -12-11j \\
-j & -1 & j & 1 \\
-1 & j & 1 & -j \\
j & 1 & -j & -1 \end{array} \right)$$

Note that this differs from Paul R's answer because of the normalization and sign conventions nobody agrees with... Divide by 4 and take the conjugate and you'll find the same result as his.


IQ imbalance amplitude I have to fix the amplitude of measured IQ receiver data. According to my information about the data, I can estimate that one signal needs to be corrected to the other in amplitude. The phase offset/imbalance can be ignored in this case. I would like to implement this in MATLAB. At the moment my I/Q signal looks like this:

I would like to correct it, that is shows a circle (would be the ideal result). But until now, I found no direct way to do that. I have no possibility to train or measure the amplitude imbalance from the I/Q receiver direclty. I have only the data. What have I to do, to normalize the Q component that the chart shows a circle?

I have already tried to calculate the amplitude offset: o=sqrt(Q^2/(1-I^2)), but this doesn't work. Both signals have a small DC offset (really small, about $10^{-4}$ in relation to values between $0$ and $1$), but this is not relevant.

What it looks like you are experiencing is carrier offset and not IQ imbalance, at least from this view that you show. Looking at the scale of the two axis it looks balanced and it is just that each axis is drawn at a different scale?

See my example plot below of QPSK before and after proper carrier recovery. Note that it also helps to show the correct symbol locations (from your timing recovery which can be done under carrier offset conditions), which are indicated by the red dots in the plot.

For more details on implementations to remove the carrier offset please see my post here:

High modulation index PSK - carrier recovery

I am interested in what particular modulation you are actually using as there may be more direct methods we could suggest, even a cross correlation of your signal with a delayed version of itself can reveal the frequency offset for example, such as depicted in the figure below), which I would be tempted to do in a post processing analysis for quickly accessing what the carrier is, that or just discerning the rate of rotation from the trajectory:

Once you successfully "stop it from spinning", here are practical techniques to remove offset errors (amp and phase imbalance) if you indeed have that problem:

First, this is what IQ imbalance would look like (in addition you could also have a "DC offset" which would appear as a stronger carrier in your spectrum):

And here is a simple technique for IQ imbalance correction:

With the techniques to establish the $\alpha$ and $\beta$ coefficients in the plots below:


Meaning of 'angle' in angle modulation I know what is angle modulation. I want to know why is named as 'Angle'. Is it because in the vector plot of sidebands are described by various angles? or what?

Well, the Angle modulation basically means adding a phase function to the carrier that depends on time.

Let's call that instantaneous phase deviation function as:

$$
\Phi(t)
$$

and our modulated signal is

$$
m(t) = A(t) * cos(\omega t + \Phi(t)).
$$

As you can see, if we eliminate the amplitude modulation (constant A(t)), we have a phase shift that varies with time in the modulated signal formula. And hey, why the heck we don't use this phase shift to transmit our information? It can be noise resistant and easy to demodulate at the same time! It's a terrific idea!

So, if we come to your question, the answer is obvious. We are modulating our information to carrier's instantaneous phase function, therefore it's angle.

If you are looking for further information and modulation techniques, you can take a look at frequency shift keying or phase shift keying modulations. Angle modulation is not limited to these modulation techniques, however they are easy to understand and give you a deeper understanding about angle modulation.

As always, have a nice day.


Is there a difference what measurement units use in covariance matrix The R matrix in the Kalman filter contains measurement noise. Diagonal elements of the matrix is the power of standard deviation. Is there a difference what measurement unit to use for standard error representation? For instance is there difference between meters and kilometers or radians and milliradians? Should all elements have a similar scale?

Most certainly yes.  The units of the states are reflected in the state covariance, and state error covariance.

The units of the measurements propagate throughout the filtering and prediction steps.

Should units be appropriately scaled? Numerical accuracy typically benefits from appropriate scaling


Weird behaviour converting velocity to displacement I'm trying to convert accelerometer data from accelerations to displacement,

I assume that the initial velocity is zero and the initial position is also zero.

Here's my code,

Fs = 500;
T = 1/Fs;
t = (0:length(y) -1)*T;

vel = cumtrapz(t, acc);
disp = cumtrapz(t, vel);

here are the respective plots from acceleration to displacement,

what am I doing wrong?

I've implemented a kalman filter but the result is the same. I'm trying to work around with wavelet decomposition, any help would be very appreciated to resolve this problem.

You have two problems here.

First, as hotpaw2 mentioned, double integration is sensitive to any non-zero bias in your accelerometer, which will certainly be present if you are using a MEMS accelerometer. Suppose you have a bias of $b$; then, your position at time $t$ is

$$p(t) = \int_0^t \int_0^{t^\prime} b \ dt^{\prime\prime} \  dt^\prime = \frac{b}{2} t^2$$.

So, your position is growing quadratically with time.  Not good!

Another problem will be the measurement noise that is riding on your signal.  Double integration of this noise signal will behave like a second-order random walk and will also cause your position to diverge over time.

A common approach to this sort of problem is to use a Kalman Filter.  Also, this is a good introduction to position estimation using accelerometers.


What is the mean power of a complex random variable? Say $\alpha$ is a complex random variable, then which one of the following expressions is correct?

$\mathbb{E}[\alpha^2]$ or
$\mathbb{E}[\alpha \alpha^*]$?

The correct expression for the mean power of a complex random variable $\alpha=x+jy$ is
$$
\begin{align}
\bar P &= \operatorname E\left[\alpha \alpha^*\right]\\
&= \operatorname E\left[x^2 + y^2\right]\\
&= \operatorname E\left[x^2\right] + \operatorname E\left[y^2\right]\\
&= \bar P_\mathrm{x} + \bar P_\mathrm{y}
\end{align}
$$

In other words, the mean power of a complex random variable is the sum of the mean powers of its real and imaginary part, respectively. In contrast, expression no. 1 from your question evaluates to
$$
\begin{align}
\operatorname E\left[\alpha^2\right] &= \operatorname E\left[x^2 + j2xy - y^2\right]\\
&= \operatorname E\left[x^2\right] - \operatorname E\left[y^2\right] + \operatorname E\left[j2xy\right]\\
&= \bar P_\mathrm{x} - \bar P_\mathrm{y} + j2\operatorname E\left[xy\right]
\end{align}
$$
Note that for the special case of $x$ and $y$ being uncorrelated $\operatorname E\left[xy\right]=E\left[x\right]E\left[y\right]$. If, in addition, $x$ or $y$ has zero mean $E\left[xy\right]=0$.


Calculating a NRR of an IIR I'd like to calculate the Noise Reduction Ratio of a system.

The system I have is $$ H(z)= \frac{z^{-5}}{1-0.1z^-1}.$$
For FIR Filters this can easily done by $$ \sum h[n].$$

For an IIR Filter this gets more complicated. I thought about making it a recursive formula for h(n) and then calculate some sort of sum which might converge (without any result yet).

Another way would be integrating $ |H(w)|$ by replacing z with $e^{jw}$. But there has to be an easier way. Some sort of trick.

I can imagine two things:

The system can be splited into two or more known systems. The NRR can than somehow be added togehter. (How would one do that?)
$z^{-5}$ is just a delay which shouldn't contribute to the NRR. The remaining part has maybe a well known NRR.

EDIT: UPDATE. I made a mistake

NRR for an FIR Filter is given by: $$NRR = \frac{\sigma^2_{Yv}}{\sigma^2_{Xv}} = \frac{1}{2\pi}\int_{-\pi}^{\pi}|H(\omega)|^2d\omega = \sum_n h(n)^2$$

EDIT 2: What's an NRR

Since I left my textbook in the office I can only show you some definitions from powerpoint slides:

  The noise reduction ratio NRR measures, how much the noise power at the output is less than the noise power at the input
  –
  The last equal sign follows from Parseval’stheorem
  –
  Note that small NRR means large noise reduction
  –
  Small NRR are good NRR

and

  $$\frac{SNR_{out}}{SNR_{in}} = \frac{E[y_s(n)^2]}{E[y_v(n)^2]} \cdot \frac{E[v(n)^2]}{E[s(n)^2]} = \frac{1}{NRR} \cdot \frac{E[y_s(n)^2]}{E[s(n)^2]}$$
  If $x_S(n)$ isn't changed by Filter $H(\omega)$ then:
  $$ \frac{SNR_{out}}{SNR_{in}} = \frac{1}{NRR} $$

Let's find the impulse response for this system:

$$H(z)= \frac{z^{-5}}{1-0.1z^-1} = z^{-5}\frac{1}{1-0.1z^-1}$$

You didn't specify the ROC, so I'll assume that we take it so that the filter is stable. Namely:

$$|z|>\frac{1}{10}$$

I've split the expression of $H(z)$ to notice that the impulse response basically consists of an exponential with a delay. We know that

$$\frac{1}{1-0.1z^-1}\xrightarrow{\mathcal{Z}^{-1}}(0.1)^nu[n]$$

where $u[n]$ is the unitary step. The $z^{-5}$ factor translates to a delay in the time domain. So the impulse response is:

$$h[n]=(0.1)^{n-5} u[n-5]$$

If we want to find the sum of all the terms squared in the sequence:

$$\begin{align}
\sum_{n=-\infty}^{\infty}h^2[n]
&=\sum_{n=-\infty}^{\infty}\left[(0.1)^{n-5} u[n-5]\right]^2\\
&=\sum_{n=5}^{\infty}\left[(0.1)^{n-5}\right]^2\\
&=\sum_{m=0}^{\infty}(0.1)^{2m}\\
&=\sum_{m=0}^{\infty}(0.01)^{m}\\
&=\frac{1}{1-0.01}\\
&=\frac{100}{99}\\
\end{align}$$

Where I've used the formula for the infinite geometric sum.

EDIT:

Note that a smaller NRR implies a worse performance of the filter at reducing noise. In this case, NRR turned out to be almost $1$ (or $0 \ \mathrm{dB}$), so the variance of the noise at the output is practically the same as in the input. As Fat32 pointed out in the comments, this depends on the position of the pole of your filter. The noise is white (as far as I understood from your comments), so it has infinite bandwidth. If the filter's bandwidth is too large, then "a lot" of noise will pass through it, appearing at the output. On the contrary, if it is narrow, then it will "filter out" most of it, leading to a decrease in the variance at the output. In this case, the filter's bandwidth is rather wide due to the pole being near the origin, and then the NRR approaches $1$: not a really good filter to reduce noise.


What is the deffinition of convexity defect in image processing? I am seeing this term in every image processing object recognition tutorial, but no one explains what this is. They only explain how to find convexity defects but what is a convexity defect?

Thank you.

In simple words, convexity defect is a cavity in an object (blob, contour) segmented out from an image. That means an area that do not belong to the object but located inside of its outer boundary -convex hull (read documentation of cv::convexHull function for details).

The image below shows it better than million of words. Areas between fingers (all marked with arrows) in this schematics of a hand contour are the convexity defects.

You can find the defects using OpenCV's function cv::convexityDefects.


Basic question regarding third octave filters I would like to use third octave filters on a time signal.
Do I have to apply all third octave filters to a single frequency or do I only need to apply the one that has its center frequency the closest to the frequency that I am looking at?

I would like to use a simple example to make my quesion clear:
Let's say I have a basic sin signal: signal = sin(2*pi*f*t) with f = 57Hz.
57Hz is almost right between the center frequencies of 50Hz and 63Hz, but the closest would be 63Hz. Do I apply all filters or only the one with a center frequency of 63Hz?

I am pretty sure that the answer is that I will always have to apply all filters in theory. Of course filters with much higher center frequencies will result in insignificant amplitudes...

Someone released my natural inheritences and regulatory systems calculated as it formed my employment, rewards  and inheritences payable to me obly all as one from more than just my g.c., Tricon, Albertsons, Tetra,  and more...as someone or a group if someones have changed my statuses to a protective state and have done so as if court ordered. I gave nevet in my life been so isolated during major life upsets such as missing children and death with zero compliance of the government and of others involved. My personally named for claim and my personal poatal mail and hijacked emails go uncorresponded and still remail displaced like my paymemt and documentations I have  applied for in the last entire seven years and my lifes worth is wromgly dispensed. From my insurance policies to my ach as named, global this state is acting as if I am dead. Even as the shut down goes, these are being processed illegally like my procured custody reformed payments from my base assigned living EXPOSURE money.
A step father and father with uncle's and aunts steal. Like they always have.
I am Lindsey ann maxfields
717 briarwood lane Hurst rd.


Denoising Oscilloscope data the right way [python3] I need to plot large oscilloscope data sets. They are too large (100 million Values) to plot all of them in a for loop. If I try to do that it results in a MemoryError. So with the help of Stackoverflow I came up with a solution that isn't pretty but it works.

Original Dataset:

To reduce the amount of data to be plotted I can remove all rows of the dataset, that have a y-value between 1 and -1. This can be done comfortably with pandas:

df is the dataset in form of a pandas.DataFrame().

df[1] is the column which is to be examined.

df.drop(df[(df[1] < 1) & (df[1] > -1)].index, inplace = True)

This results in:

This reduces the length of the dataset from 100 million to about 200 thousand. My Question is: Is there a way to reduce data length and keep its original look (prevent the ugly cut-off)?

This stack seems to work differently than Stackoverflow. Pardon me for all wrongdoings.

If you are simply interested in plotting the data then any data reduction technique would do,even if it appears to be crude.

Effectively, the plotting function itself will not plot all the data, because the space assigned to the plot has a finite number of $N_x \times N_y$ pixels assigned to it. For example, if your plotting area was $1024 \times 768$, then any waveform with more than 1024 samples would need to somehow be re-scaled or re-sampled. Plotting functions do that internally anyway but always retain a link to the original data to do the whole process again, usually onUpdate.

What you could try to do is to resample your waveform and specifically downsample it.

Downsampling is composed of two steps: The step of anti-alias filtering and the step of decimation.

If you are purely interested in plotting a waveform, you can try to simply decimate it with a low factor (say, retaining 1 sample out of every 5) and then send it for plotting. It will still have more data than required by the plot but it will be of more manageable size.

The reason for opting for pure decimation here is because the filtering step (although necessary) can be costly (in terms of computation) and take a bit longer. If you can afford the filtering time however, I would suggest that you do a proper downsampling.

If you are using Python, you can try this function (or this one) to handle the downsampling step.

It goes without saying here that if you wish to process the signal further, you must downsample it properly.

Finally, another reason why you are running out of memory might be the fact that you are retaining all data in memory even when you are done using it. Therefore, once you have sent a waveform for plotting, try del()eting it, possibly followed by gc.collect() to ensure that that part of memory will become immediately available.

Hope this helps.


Calculating SNR & Es/No and placing noise onto a signal I have a signal power of -12 dBW from a 16 QAM signal after the first RRC.
I want to create a channel SNR of 30 dB
I calculate my noise power from this, i.e. its -42 dBW and make the noise to add to signal
#Power values below have been converted to Watts
noise_distribution = np.random.randn(N)  / np.sqrt(2)
noise = np.sqrt((noise_power/signal_power)) * noise_distribution
signal = signal + noise

This noise is inserted into the signal after the first RRC and before the second.
I know my signal power is spread over a symbol rate of 2 MHz, with a 0.35 Roll off on RRC.
Signal PSD = -12 - 10log10(2MHz(1+(2*0.35))) = -77 dBW/Hz
However my noise power is spread over Fs:
-32 dBW - 10log(40MHz) = -95 dBW/Hz
Is this the correct way to do SNR or have I just done Es/No because it doesn't seem to work in my mind, it doesn't end up being 30 dB in PSD terms.

I assume the OP wants the SNR to be 20 dB just prior to decision in the receiver, in which case the noise would be 20 dB below the signal power AFTER being passed through the second RRC filter. The total noise in the bandwidth of this filter would be 20 dB below the total noise of the signal. Since neither if flat over bandwidth, I would not recommend units of dBW/Hz, but work with total power (or as the square root; the standard deviation of signal and noise).


How to understand the implementation codes in Debevec and Malik's HDR paper The paper titled "Reconvering High Dynamic Range Radiance Maps from Photographs" from Debevec and Malik is a seminal paper in high dynamic imaging. In this paper, it proposed a method to estimate the camera response function (CRF), which is equal to optimize the following objective function:

$$\mathcal O = \sum_{i=1}^N\sum_{j=1}^P\left\{w\left(Z_{ij}\right)\left[g\left(Z_{ij}\right)-\ln E_i-\ln\Delta t_j\right]\right\}^2+\lambda \sum_{z=Z_{\rm min}+1}^{Z_{\rm max}-1}\left[w(z)g''(z)\right]^2$$

The author also provides MATLAB codes to optimize this object function:

%
% gsolve.m − Solve for imaging system response function
%
% Given a set of pixel values observed for several pixels in several
% images with different exposure times, this function returns the
% imaging system’s response function g as well as the log film irradiance
% values for the observed pixels.
%
% Assumes:
%
% Zmin = 0
% Zmax = 255
%
% Arguments:
%
% Z(i,j) is the pixel values of pixel location number i in image j
% B(j) is the log delta t, or log shutter speed, for image j
% l is lamdba, the constant that determines the amount of smoothness
% w(z) is the weighting function value for pixel value z
%
% Returns:
%
% g(z) is the log exposure corresponding to pixel value z
% lE(i) is the log film irradiance at pixel location i
%
function [g,lE]=gsolve(Z,B,l,w)

n = 256;

A = zeros(size(Z,1)*size(Z,2)+n+1,n+size(Z,1));
b = zeros(size(A,1),1);

%% Include the data−fitting equations

k = 1;
for i=1:size(Z,1)
  for j=1:size(Z,2)
    wij = w(Z(i,j)+1);
    A(k,Z(i,j)+1) = wij; A(k,n+i) = −wij; b(k,1) = wij * B(i,j);
    k=k+1;
  end
end

%% Fix the curve by setting its middle value to 0

A(k,129) = 1;
k=k+1;

%% Include the smoothness equations

for i=1:n−2
  A(k,i)=l*w(i+1); A(k,i+1)=−2*l*w(i+1); A(k,i+2)=l*w(i+1);
  k=k+1;
end

%% Solve the system using SVD

x = A\b;

g = x(1:n);
lE = x(n+1:size(x,1));

The implementation codes are well understood, and my question is why this implementation can optimize the object function. If I understand well about this implementation, the solution to the object function is the one that makes the gradients of the object function equal to zero. However, one thing I am confused is that why the first part and the second part of the object function are treated independently.  How could it be possible?

OK, this is incomplete but this is how far I've gotten:

My previous answer, which you can see in the edits was incorrect.

Your question boils down to:

  However, one thing I am confused is that why the first part and the second part of the object function are treated independently.

So let's try having a look at the objective function:
$$
\begin{array}{rcl}
{\cal O} &=& \displaystyle \sum_{i=1}^{N} \sum_{j=1}^P \left\{
w(Z_{ij}) [ g(Z_{ij}) - \ln E_i - \ln \Delta t_j ]
\right\}^2 + \lambda \sum_{z =Z_{\rm min}+1}^{Z_{\rm max}} [w(z) g^"(z)]^2\\
&=& \displaystyle \sum_{i=1}^{N} \sum_{j=1}^P w(Z_{ij})^2
 [ g(Z_{ij})^2 + (\ln E_i)^2 + (\ln \Delta t_j)^2
- 2 g(Z_{ij}) \ln E_i - 2 g(Z_{ij}) \ln \Delta t_j + 2 \ln E_i \ln \Delta t_j
 ]\\
&&+ \lambda \displaystyle \sum_{z =Z_{\rm min}+1}^{Z_{\rm max}} [w(z) \{ g(z-1) - 2g(z) + g(z + 1) \} ]^2
\end{array}
$$
where the second line just expands the first term and substitutes for $g^"$ as the paper suggests.

Now, we need to choose the $\ln E_i$ and $g$ values to minimize this so we must take the derivative of $\cal O$ with respect to $\ln E_i$ and to $g$.

$$
\frac{d \cal O}{d \ln E_i} = \sum_{i=1}^{N} \sum_{j=1}^P w(Z_{ij})^2
 [  2 \ln E_i
- 2 g(Z_{ij}) + 2 \Delta t_j
 ]
$$