Is applying a linear operation to the input and output of a continuous linear operator OK? The scientist and engineers guide to digital signal processing, Steven W. Smith, p. 134-135 states:

  [...] imagine a linear system receiving an input signal, $x[n]$. and generating an output signal, $y[n]$.  Now suppose the input signal is changed in some linear way, resulting in a new input signal, which we will call $x'[n]$.  This results in a new output signal, $y'[n]$.  The question is, how does the change in the input signal relate to the change in the output signal?  The answer is: the output signal is changed in exactly the same linear way that the input signal was changed. [...] A linear change made to the input signal results in the same linear change to the output signal.'

I assume this can be written for the continuous case as:

If $g = L(f)$, then $M(g) = L(M(f))$

where $L, M$ are linear (LTI) operators, $f, g$ are functions corresponding to $x[n], y[n]$ and  $M(f), M(g)$ correspond to $x'[n], y'[n]$.

My question is: is this true for the continuous case?

Any references, or a proof would be appreciated.

This question is a near duplicate of this question (answered) and also this one (now deleted).  I am sorry about this-spent much time on the web but only found the reference I quoted-but really need an answer asap.

To be a bit picky, the claim does not hold for all linear systems. Here's a counter-example:

Let the system $Tr$ be defined by the input-output relation $$y[n]=Tr\{x[n]\}=nx[n].$$
This system is clearly linear, since

$$Tr\{ax_1[n]+bx_2[n]\}=n(ax_1[n]+bx_2[n])=aTr\{x_1[n]\}+bTr\{x_2[n]\}$$

Now, let's apply the following linear operation $H$ to an input signal: $y[n]=H\{x[n]\}=x[n-1]$. This operation is also linear. Let's see what happens to the claim:

$$Tr\{H\{x[n]\}\}=nx[n-1]\neq H\{Tr\{x[n]\}=(n-1)x[n-1].$$

So, in general the claim does not hold for a linear system.

However: Assuming the claim is stated for linear time-invariant systems (LTI), then it is correct. (The claim was missing the constraint of time-invariance). The proof is simple: Any LTI system is described by its impulse response, and the output of a system is the convolution of the input signal with its impulse response. Given a system $L$ with impulse response $l[n]$, and a system $M$ with impulse response $m[n]$, we get

$$g=L\{f\}=l*f$$
$$M\{g\}=m*g = L\{M\{f\}\}=l*m*f=m*l*f=M\{L\{f\}\}$$

The reason why this is true is the commutativity of the convolution operation (i.e. $x*y=y*x$). Since the commutativity holds for both discrete and continuous time, the claim holds for both discrete and continuous time.


Calculating impulse response of system Given is the illustrated circuit diagram of a linear, time-invariant, time-discrete system:

How do I show that the total system has the impulse response $h[n] = aδ[n] + bδ [n - 1] + cδ [n - 2]$ and determine the constants a,b and c?

with the following impulse responses of the subsystems:

$h_1[n] = δ[n] + (1/2)δ[n − 1]$

$h_2[n] = δ[n] − (1/2)δ[n − 1]$

$h_3[n] = −δ[n]$

$h4[n] = − (1/2)^n σ[n]$

So what I have done:

$h_{3||4}=−δ[n]−(1/2)^n σ[n]$

and than:

$h_{3||4+2}=\sum_{k=-\infty}^{\infty}(δ[n-k] − (1/2)δ[n-k − 1])(−δ[n]−(1/2)^n σ[n])$

$h_{3||4+2}=−δ[n]−(1/2)^n σ[n]-\frac{1}{2}(−δ[n-1]−(1/2)^{n-1} σ[n-1])=−δ[n]-\frac{1}{2}^nδ[n]+\frac{1}{2}δ[n-1]$

$h_{3||4+2+1}=(\frac{1}{2})^nδ[n]-δ[n-1]$

But the solution is $a=1$ and $b=-1$. How is $a=1$?


Should I consider rgb values of a pixel as one value? I was reading this paper for a project work.

  We train on 1.6 million 32*32 color images that have been preprocessed
  by subtracting from each pixel its mean value over all images and then
  dividing by the standard deviation of all pixels over all images.

I've trouble distinguishing between "from each pixel its mean value over all images" and "standard deviation of all pixels over all images".

Since, I'm dealing with color images, can I just take rgb values of each pixel as one value or should I calculate the mean and SD for every color differently?

For example if I have r=255, g=255, b=255, can I take pixel value as (in binary), (r<<16)+(g<<8)+b ?

For example if I have r=255, g=255, b=255, can I take pixel value as (in binary), (r<<16)+(g<<8)+b ?

While this number would identify the colour, it would not help with a mean or standard deviation. e.g. pure blue would be 255 but pure red would be 16711680

The paper is unclear so you have a few choices:

Process each channel individually
Convert the image to grayscale, find the mean and std dev, and use this for each channel.

There may be other methods that use different colour spaces.


Phase is not linear in frequency domain after delay in time domain When impulse response is symmetric around sample zero phase entirely zero.

in = 0 , 0 , 0 , 0 , 1 , 1 , 1 , 1 , 0 , 0 , 0 -> This gives zero phase result

As I know if I apply a delay in time domain this should result a linear phase in frequency domain. So if I delay my previous vector by one;

in = 0 , 0 , 0 , 0 , 0 , 1 , 1 , 1 , 1 , 0 , 0 -> This should give linear phase responce

const int N = 11;

// in = 0 , 0 , 0 , 0 , 0 , 1 , 1 , 1 , 1 , 0 , 0
std::vector< std::complex<double> > in (N);
std::vector< std::complex<double> > out (N);
std::vector< std::complex<double> > polarOut (N);

auto middleElem = in.begin() + in.size()/2;
//std::fill( middleElem - 1, middleElem + 3, 1); //zero phase
std::fill( middleElem, middleElem + 4, 1);

fftw_plan my_plan  = fftw_plan_dft_1d(N, reinterpret_cast<fftw_complex*>(&in[0]),
                                         reinterpret_cast<fftw_complex*>(&out[0]), FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(my_plan);

std::transform( out.begin(), out.end(), polarOut.begin(),
                []( auto& in ){
                                return  std::complex<double>( std::abs(in), std::arg(in) );
                              } );

The result of the fourier transform in polar notation is :

    polarOut    <11 items>  std::vector<std::complex<double>>
        [0] (4.000000, 0.000000)    std::complex<double>
        [1] (3.228707, 2.570394)    std::complex<double>
        [2] (1.397877, -1.142397)   std::complex<double>
        [3] (0.372786, -1.713596)   std::complex<double>
        [4] (1.088156, 0.856798)    std::complex<double>
        [5] (0.546200, -2.855993)   std::complex<double>
        [6] (0.546200, 2.855993)    std::complex<double>
        [7] (1.088156, -0.856798)   std::complex<double>
        [8] (0.372786, 1.713596)    std::complex<double>
        [9] (1.397877, 1.142397)    std::complex<double>
        [10]    (3.228707, -2.570394)   std::complex<double>

What am I missing here? Why phase responce is not linearly increasing?

The phase is in fact linear, apart from jumps of $\pi$ at the zeros of the magnitude. There are two reasons why you don't see it:

the phase you compute is the principal value, which is in $(-\pi,\pi]$. Note that you can always add or subtract a multiple of $2\pi$ to the phase without changing anything. Removing these artificial jumps due to the principal value is called phase unwrapping.
you only compute very few points, and due to point 1 above you just see an erratic behavior.

Solution:

compute a more detailed DFT (by zero-padding the original signal)
unwrap the phase (if you have Matlab/Octave, there's a command to do that)

Here an example in Matlab/Octave:

x = [0 0 0 0 0 1 1 1 1 0 0];
X = fft(x,512);         % length 512 FFT
ph = arg(X);            % wrapped phase
phu = unwrap(ph);       % unwrapped phase
plot((0:256),phu(1:257)/pi), grid on


Why does this plot correspond to this function? I have a function:

$$
u_1(\mathbf{x})=
\begin{cases}
1, & \text{if} \;\;|x_1|\le r_1, |x_2| \le r_2 \\
0, & \text{otherwise}
\end{cases}
,\mathbf{x}=[x_1,x_2]^T \in \mathbb{R}^2
$$

And I'm told that if $r_{1}=r_{2}=5$, then plotting it in Matlab with stem3 should look like this:

Maybe I'm just too exhausted at this point, but I don't see how this is the same function. It looks like it's equal to 1 between 10 and 20 rather than between -5 and 5? Can someone please explain what's going on here?

The axes in your plot seem to be indexes into an array, as if you were using stem3(Z) instead of stem3(X,Y,Z).

Check the example in the documentation.


Wavelet transform of a spatial convolution Does anyone know if there exist a kind of convolution theorem for the discrete wavelet transform (decimated or undecimated)?

In other words can I find a simple form of
$W\left[ \int f(t) g(x-t) \, dt\right] $ where $W$ is the discrete wavelet transform operator?

I cannot say I have a clear understanding of this at this time.  However, a few pointers. I'd love to see somebody provide a detailed account. Others bits at: Multiplication in the wavelet domain, what does it look like in real space?

The nonexistence of a wavelet function admitting a wavelet convoluton theorem of the Fourier type, 1994, A. R. Lindsey, unpublished report
Convolution using the undecimated discrete wavelet transform, 1996

  Convolution is one of the most widely used digital signal processing
  operations. It can be implemented using the fast Fourier transform
  (FFT), with a computationalcomplexity of $O(N \log N)$. The
  undecimated discrete wavelet transform (UDWT) is linear and shift
  invariant, so it can also be used to implement convolution. In this
  paper, we propose a scheme to implement the convolution using the
  UDWT, and study its advantages and limitations.

Convolution theorems for linear transforms, 1998

  This correspondence explores the existence of convolution theorem for
  linear transformations under a variety of different assumptions. There
  are eight convolution theorems, all Fourier-related with only N
  operations in the transform domain and no ordering constraints on the
  convolution components in the result. They include circular
  convolutions and correlations.

The convolution theorem for the continuous wavelet tranform, 2004

  We study the application of the continuous wavelet transform to
  perform signal filtering processes. We first show that the convolution
  and correlation of two wavelet functions satisfy the required
  admissibility and regularity conditions. By using these new wavelet
  functions to analyze both convolutions and correlations, respectively,
  we derive convolution and correlation theorems for the continuous
  wavelet transform and show them to be similar to that of other joint
  spatial/spatial–frequency or time/frequency representations. We then
  investigate the effect of multiplying the continuous wavelet transform
  of a given signal by a related transfer function and show how to
  perform spatially variant filtering operations in the wavelet domain.
  Finally, we present numerical examples showing the usefulness of
  applying the convolution theorem for the continuous wavelet transform
  to perform signal restoration in the presence of additive noise.

On the uniqueness of the convolution theorem for the Fourier transform, 2008

  This paper shows that members of the fourier transform family are the
  only linear transforms that have a convolution theorem, that is, that
  can replace $O(N^2)$ operations of a convolution in a time domain by
  $O(N)$ operations in a transform domain. Generally, there is an
  additional cost to compute the transform itself. Our observation is
  motivated by recent activity in wavelet and subband decompositions and
  related spectral analyses, which are attractive alternatives for
  signal compression applications. A natural question when using such
  techniques is to determine if convolutions of $N$-point signals can be
  calculated with fewer operations in a compressed transform domain than
  in an uncompressed time domain. The answer is negative for a broad set
  of assumptions. This paper indicates what assumptions must be relaxed
  in seeking a linear transform that has a convolution theorem
  comparable to the convolution theorem for Fourier transforms.


Do I get less accurate results if I do not use a power of 2 points for my FFT? I heard that the padding with zeros is mainly for efficiency and speed, other than that, is there a downside to use fft(signal) instead of fft(signal, N) where N is a power of 2 assuming signal length is not a power of 2.
Thank you

Assuming a proper implementation, and discarding the small rounding differences that you'll get just by doing a different sequence of operations, no, there's no substantive drop in accuracy if you slightly pad your DFT length to the next power of 2. Whether it ends up actually being any faster is dependent upon your implementation, however. Why is that?

It's a very common misconception that there is a single Fast Fourier Transform algorithm and that you must use a power-of-2 size for best performance. If you look for a description of an FFT algorithm, you'll often see the radix-2 decimation-in-time or decimation-in-frequency techniques explained, likely because they're the easiest to illustrate. However, even the seminal FFT technique, the Cooley-Tukey algorithm, generically factorizes the FFT size into smaller numbers, not just powers of 2.

Using a good FFT library, you'll get the best performance if your FFT size can be factored into a number of small prime factors. The FFT library will then have optimized implementations of DFT kernels for each of these primes, which can then be recombined appropriately to yield the full set of DFT outputs. As I alluded to in a comment on another question, modern libraries will often give good performance for all prime factors ~13 and below.

With that said, you may find that by padding up to the next power of 2, it's possible that your transform might become slower. If you're already using an FFT size that your library implementation is well-suited for, you're just adding extra work for yourself by padding the size out. The best way to judge this is to benchmark a few candidate sizes and see which does best on your platform. If you have to make an automated choice of a good FFT size, then based on the characteristics of the radixes that your library supports, you can choose the next size that has an appropriate set of prime factors.


Sinogram Fourier Analysis Doesn't Show Clear Bow Tie Shape I'm working on finding the center of rotation of a set of test tomographic projections in order to perform 2D reconstruction. I'm trying to implement the  algorithm to find center of rotation as described in Vo et al. 2014, which involves Fourier transforming the sinogram and searching frequency space outside the "bow tie" double-wedge region for minimum artifacts that result from wrong center of rotation. However, when I Fourier transformed a sample sinogram generated from Shepp-Logan model, the resulting frequency space has lots of artifacts outside the double-wedge region even when center of rotation is correct. Through visual inspection I've also noticed that the amount of artifacts are sensitive to number of projections in the sinogram.

Below: FT of sinogram with 200 projections, from 0 degree to 180 degree

Below: FT of sinogram with 180 projections, from 0 degree to 180 degree

The FFT'ed sinogram in the paper shows distinct zero regions outside the double wedge area, which is mathematically true as it is proven here. However, I'm not sure why my FFT'ed sinograms do not show similar zero regions as in the paper.

This is my first time posting here, so please let me know what other information is needed for others to help me. Thank you.

Below is my code

P=phantom(256);

angs=linspace(0,179,200);

pjs=radon(P,angs);

pjs=rot90(pjs,3);

pjs_inv=fliplr(pjs);

pjshift_inv=circshift(pjs_inv,[0 0]);

sinogram=[pjs; pjshift_inv];

figure(1),imagesc(sinogram);axis image,colormap gray

ft=real(fftshift(fft2(sinogram)));

figure(2), imshow(ft),colormap gray,axis image,axis on

Try imtool() and adjust the contrast to show ft similar as in the paper. I set eliminate outliers to 5%, then (contrast) window minimum to 0.0. May I suggest to ask the author about the settings for the display range in his paper.

Now to the signal processing.

The mentioned paper by Vo is public available. Moreover Vo published the algorithm in mathematica code. Note that the code by Vo does a gaussian filtering, i.e. a lowpass of the sinogram. IMHO to achieve similar results as in a real scan.

Generally a real scan shows lowpass behaviour, i.e. the high frequencies in the projections in the sinogram are lost due to the size of the focal point inside the x-ray tube, the detector system and radiation scatter inside the object. In contrast the forward projection algoritm: Those may keep high frequencies in the projections, depending on their base algorithm.

sinogram1 = GaussianFilter[ sinogram1, {{3,2}}];
(* Apply smooth filter: recommend blur more in the vertical direction *)

The Mathematica help tells, that those 3 and 2 steer the xy-sigma of the GaussianFilter(). One can go theroretically through that to find the values of the matrix, but I prefer to run the above filter over a dirac impulse to see the resulting filtermatrix and establish it similar in matlab for use in imfilter(). But I have no access to mathematica. May be s.o can translate the Mathematica GaussianFilter() into a matlab fspecial().

So let's filter the sinogram with a lowpass to ged rid of the high frequencies introduced by radon() :

imtool(ft); %% eliminate 5% outliers and set min to 0.0
h=fspecial('gaussian',[2 3]); %% wider in projection direction
sino_lp = imfilter(sinogram,h,'replicate');
ft_lp=real(fftshift(fft2(sino_lp)));
imtool(ft_lp);

Such a fourier transformed sinogram has lower values in the high frequencies.


Why is the received signal=information+RF signal+noise,not just information+noise? In this paper: https://www.researchgate.net/publication/325008791_Beamforming_Design_for_Max-Min_Fair_SWIPT_in_Green_Cloud-RAN_with_Wireless_Fronthaul

Here is its system model

In this model,it said the received signal in each DR is

$y_k=\mathbf h_k^H \mathbf w_k s_k^D+\sum\limits_{i=1,i \neq k}^{K}\mathbf h_k^H \mathbf w_i s_i^D+\sum\limits_{j=1}^{J}\mathbf h_k^H \mathbf v_j s_j^E+n_k$, Notice that due to the energy symbol $s^E_j$ is randomly generated, which carries no information but only satisﬁes the RF regulations;$s_k^D$ is data symbol for DR k .

Honestly,i don't understand about $s^E_j$,i mean i know the received signal of DR1,

$y_{DR1}=$(signal from RRH1) + (signal from RRH2 and RRH3,etc)+(noise),that is

$y_k=\mathbf h_k^H \mathbf w_k s_k^D+\sum\limits_{i=1,i \neq k}^{K}\mathbf h_k^H \mathbf w_i s_i^D+n_k$

but i don't know where does  $\sum\limits_{j=1}^{J}\mathbf h_k^H \mathbf v_j s_j^E$ come from

I have a thinking before,the signal RRH transmit is information+RF signal,that is $s_k^D+s_j^E$,but in the other paper,they will only mention the $s_k^D$,like this paper: https://arxiv.org/pdf/1805.08898.pdf ,it just said $y_k=\mathbf h_k^H \sum\limits_{j=1}^K \mathbf f_j s_j+n_{a_k}$,so does anyone know about what the $s_j^E$ is?

Well, it seems the authors assume a setup where there are separate data receivers (DR) and energy receivers (ER) ("To achieve a balanced user experience for separately located datareceivers (DRs) and energy receivers (ERs) in the network, joint transmit beamforming vectors are optimized [...]"). In that setup, $s_k^E$ would be symbols that are used for sending energy to the ERs whereas $s_k^D$ are data symbols for the DRs.

Therefore, if you are a DR and want to receive data, you care about two types of interferences: data symbols for other DRs and energy symbols for the ERs. Note that you didn't copy equation (1) exactly: the second term should be $s_i^D$ (data symbol interference). Then you have the third term with interference of energy symbols for the ERs and the fourth term for the noise.

As an ER, you wouldn't care about interference, it's actually additional power, if you like.

Other papers use other models where one node may want to receive data and harvest energy at the same time. As they say: "all models are wrong, but some are useful..."


An invertible system with memory Suppose $\mathcal{L}$ be invertible system with memory. Does $\mathcal{L}^{-1}$ have memory necessarily?

Intuitively I think the answer is "yes". There are many examples showing that. For instance $\mathcal{L}(x(t)) = x(t-2)$ and $\mathcal{L}(x(t)) = x(\frac t 3)$. Another example which seems problematic to me is
$$\mathcal{L}(x(t)) = \int_{-\infty}^{t}x(\lambda)d\lambda$$The inverse is $$\mathcal{L}^{-1}(x(t)) = \frac{dx(t)}{dt}$$Does differentiator have memory? Of course the main question here is about memory of an invertible system which has memory. Note that here $\mathcal{L}$ can be nonlinear as well.
For clarity, I add some related definitions from Oppenheim's book:
Invertible system: A system is said to be invertible if distinct inputs lead to distinct outputs.
Causal system: A system is causal if the output at any time depends only on values of the input at the present time and in the past.
Memoryless system: A system is said to be memoryless if its output for each value of the independent variable at a given time is dependent only on the input at that same time.

For time-based systems, I understand that it is difficult to imagine a memory of the future. But for general systems, $-t$ and $t$ are just left and right (think of a spatial system). Other discussions are in LTI system $y(t)=x(t−T)$ with or without memory, What is a memory less system?, or A question about the concept of the time.
By definition of invertibility, $\mathcal{L}^{-1}$ is such that  $\mathcal{L}^{-1}( \mathcal{L}(x))=x$. But also that $\mathcal{L}( \mathcal{L}^{-1}(x))=x$ (by the way, derivatives and integrals are not inverses). Let us suppose the converse: $ \mathcal{L}^{-1}$ has no memory. Hence $\mathcal{L}^{-1}(x[n]))$ can only use the present state, and $\mathcal{L}$ as well to yield $(x[n])$.
So, if $\mathcal{L}^{-1}$ is memoryless, $\mathcal{L}$ is memoryless as well. By contraposition, the converve is true


Finding correct values based on information from two arrays Consider the following scenario:

Say, one machine is sending out a beep signal every 10 seconds in a very noisy environment. I have two sensors which detects these beeps independently. Device A is near to machine compared to B, so B always have a two samples delay compared to A. Since it is noisy environment, devices sometimes detects noise also as beeps and sometimes even misses the data.

Below is the data

A = [5,10,20,25,30,40,45,50]
B = [12,18,22,32,38,42,52]

See, A and B have beeps at 10 & 12 which says beep @5 in A is a noise. Similarly @25 in A because 30 in A and 32 in B says they are the good data.

Similarly in B, 18 is a noise because there is a synchronized peak at @20 in A and @22 in B.

Below is a plot of that (red mark shows noisy beeps):

Question:

Is there any mathematical/signal processing method which I can use to find the location of correct beeps based on the synchronized information available from two datas and get rid of the noisy ones?

This is a fairly easy task. The reason is that your signals are monotonically increasing and encode the same additively evolving time information.
First of all, I suggest you subtract 2 from B. Then what you do is:

1) Start at position 0 on A and position 0 on B and maintain 2 separate indices (say $i$ and $j$) for traversing A and B.

2) Compare current values of $A$ and $B$: $A_i$ and $B_j$.

3) The moment where $A_i > B_j$ increase $i$ : $i=i+1$.

4) The moment where $B_j < A_i$ increase $j$ : $j=j+1$.

5) The moment where $B_j = A_i$ you have a beep. Increase both indices.

6) Goto step 2 and repeat until you reach the end.


Detect ending of frame in H.264 video stream I am using JM reference software for encoding and decoding of video. What I want to do is to add Redundant bits using RS codes to the encoded video so that if some bits get lost during transmission we recover it on decoder side.

The problem I am getting is to detect end of frame when reading H.264 file.
So my questions are as follows

How to detect end of frame when reading h.264 file?
Also, how to detect end of GOP in h.264 stream?

I think Technically speaking the only way to detect when a frame ends is to see where the next one begins. That means, you should search for the next start prefix code 0x000001. If you go to the previous data before this (which in the general case is the byte alignment syntax elements) and add your RS codes, a compliant decoder should ignore all that until the prefix code is seen.

As for the end of GOP, I'm afraid that the GOP concept does not exist in the Video standard, in theory you could start receiving frames which use previous frames as reference, but there is never a repeating structure. You can detect the end of a coded video sequence, but that does not necessarily coincide with what we know in the video community as "the GOP".

Edit:
In reply to the comment:
Lets suppose the bit stream is this (I am going to invent, it will be all garbage):
0x 00 00 01 e5 d8 31 04 ab 80 00 00 01 8e 56 04 94 cd 80 10 00 00 01 ...

You know for example that 00 00 01 e5 d8 31 04 ab 80 is the first frame, so you take all those bits and generate a CRC for them (or any other error correcting technique that you want), lets suppose you insert one byte of CRC ( that hypothetically will allow you to detect up to 5 bits in error and correct 3 for example), then your new bit stream would be

0x 00 00 01 e5 d8 31 04 ab 80 fc 00 00 01 8e 56 04 94 cd 80 53 10 00 00 01

here the CRC would be the bytes in boldface.

As for the GOP. Remember that even though when you decode you can see a repeating structure of frames "it is not required" to be like that forever. IF you want you could count the number of frames between I frames and in the practical cases it will work. Most H.264 encoders use a structure like I  {[B] (xN) P} (xM) I ...

Examples:

I B B P B B P B B I B B P B B P B B I B B P B B P B B ... etc

I B B B P B B B P B B B I B B B P B B B P B B B I B B B P B B B P B B B ... etc

But if you start counting and see in first example that the number of frames between two I frames is 8, nothing guarantees that in the next cycle it would be the same. Usually it will be, but there is nothing in the standard saying that you have to do it like that.

Hope it helps


How is the Laplacian used for Unsharp Mask? Some unsharp algorithms (e.g. Matlab's fspecial) use a Laplacian kernel to achieve sharpening of the image.

From what I know, the image needs to be blurred first, and then the blurred image is subtracted from the original to get a sharper image. The Laplacian part comes exactly where alleged smoothing should take place. Correct me if I'm wrong, but Laplacian doesn't do any smoothing.

So how does the Laplacian help in sharpening the image?

The regular sharpening i.e. high pass filter, essentially applies the inverse filter of blurring.

Basically, when you do unsharp masking - instead of subtracting the blurred image, you add the edge image as a portion of it. we can have $g(x,y)$

$$g(x,y)=f(x,y)-f_{smooth}(x,y)$$

which is nothing but an edge image using LoG.

The sharp image is produced by adding portion of this is added back to original image and hence contrast near edge increases drastically. This is how you produce the sharp image $f_{sharp}$

$$f_{sharp}(x,y)=f(x,y)+k*g(x,y)$$

Courtesy http://homepages.inf.ed.ac.uk/rbf/HIPR2/unsharp.htm.

More details are given in that page.


How can I find expansion coefficients of the y vector in a given basis? Consider the following vectors in $\mathbb R^4$:

$$\mathbf{v}^{(0)}=\begin{bmatrix}\frac{1}{2}\\\frac{1}{2}\\\frac{1}{2}\\\frac{1}{2} \end{bmatrix} , \mathbf{v}^{(1)}=\begin{bmatrix}\frac{1}{2}\\\frac{1}{2}\\-\frac{1}{2}\\-\frac{1}{2} \end{bmatrix}, \mathbf{v}^{(2)}=\begin{bmatrix}\frac{1}{2}\\-\frac{1}{2}\\\frac{1}{2}\\-\frac{1}{2} \end{bmatrix}
\mathbf{v}^{(3)}=\begin{bmatrix}\frac{1}{2}\\-\frac{1}{2}\\\frac{-1}{2}\\\frac{1}{2} \end{bmatrix}
$$

Let $$\mathbf{y}=\begin{bmatrix}0.5\\1.5\\-0.5\\0.5\\\end{bmatrix} \text,$$

what are the expansion coefficients of $\mathbf{y}$ in the basis $\{\mathbf{v}_0, \mathbf{v}_1, \mathbf{v}_2, \mathbf{v}_3\}$?

You need to find projection of $\vec{y}$ along unit vectors in the direction of each of the basis vectors $\mathbf v^{(i)}$.

For finding unit vector in the direction of the vector, you just divide the vector by its magnitude.

And, for finding projection along a unit-vector, you just take the dot-product with the unit-vector.

So, the combined steps boil down to doing the following: $$<\vec{y}, \frac{\mathbf v^{i}}{||\mathbf v^{i}||^{2}_{2}}>$$

In your case the basis-vectors are already normalized, meaning they are already unit vectors, so just take the dot-product(inner product) of vector $\vec{y}$ with each of the basis vectors $\mathbf v^i$


Samples per second in a row vector in Matlab? I wrote two functions in MATLAB, function called half that removes every other element from an arbitrary length vector, creating a shorter vector made of only the odd-numbered elements of the original vector; and a MATLAB function double that creates a longer vector by adding an additional element between neighboring elements in the original vector (Each new
element should equal the average of its neighboring elements).
If you sample a $100\textrm{ Hz}$ tone at $8000$ samples per second, and store the samples as a row vector, what will be the effect of half and double on the pitch and duration of the tone if the playback rate is identical to the sampling rate? In other words, what is the new pitch in Hertz, and duration in seconds of the new sample in each case?

Duration is just the length of the signal. half makes the duration half, and double makes it double!
Pitch is the auditory sensation of the frequency. If you double the length of a signal (given the rate of playback remains the rate of sampling) then it will sound like a lower frequency. Similarly, when you make it half, it will sound like a higher frequency.

The code you can use to hear it:

fs = 8000;  % sampling frequency
duration = 1; % seconds
f0=100; % frequency
t=0:1/fs:duration; %time axis
wave = 10*sin(2*pi*f0*t);
sound(wave,fs) % original tone
pause
half_wave = wave(1:2:end);
sound(half_wave,fs) % half duraton, double pitch
pause
double_wave = interp1(linspace(0,1,numel(wave)),wave ,linspace(0,1,numel(wave)*2-1));
sound(double_wave,fs) % double duration, half pitch


How to implement Signal transmitted by RF source I am implementing a research paper, in which i have to implement signal transmitted by RF source in matlab. The equation given in this paper is : $s(n)e^{j \cdot 2 \cdot \pi \cdot f \cdot n}$. Any help in this regard would be highly appreciated.

I assume this would be as follows:

s*exp(j*2*pi*f*n)

Where:

n: sample index
s: baseband signal, magnitude and phase for each sample
f: carrier frequency in cycles/sample

Note how the carrier as a discrete signal is given in cycles per sample. This is the normalized frequency and is related to the desired frequency $F$ in Hz by dividing by the sampling rate $f_s$ as follows:

$f = F/f_s$


Equalizing audio under linux / Raspberry Pi I have an audio streaming receiver written in Python running on a raspberry pi. This receiver writes the streamed audio out to an alsa device, which in turn passes it to an amplifier which is driving a DML speaker that I made. DML speakers require rather heavy equalization to make them sound correct - so I want to add in digital signal processing on the rpi.
Using Room EQ Wizard, I have created a correction filter. I am able to see the raw equalizer values (frequency, gain and Q), and also export the filter as an impulse response wav. I have tested this IR wav using shairport-sync (I even wrote the directions) and it works very well.
The question is how to implement the equalization under linux (headless, no GUI) for audio that is coming from a python script.
I think there are three different routes that could work -

using a software based graphic equalizer
applying convolution to the audio signal either under ALSA after the audio is written by python
or directly in python before it is written

There is an alsa equalizer plugin, however this is too limited as there aren't enough bands, the bands are preset frequencies, and the Q is not adjustable.
There is also an experimental alsa convolution plugin, which unfortunately I can't get working.
I have also found this question which lead me to ponder whether or not python could apply the convolution directly. Unfortunately I am not a python (or audio processing) expert, so this would take quite a lot of effort to figure out.
Is anyone aware of an equalizer or alsa plugin that will do what I need?

You can simply implement this as an IIR filter directly in Python using https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.sosfilt.html
Each Room EQ wizard section can be turned into a second order filter section using RBJs Audio EQ Cookbook (See for example https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html or https://www.musicdsp.org/en/latest/Filters/197-rbj-audio-eq-cookbook.html). You want the "peaking EQ". Then just cascade all sections into a single SOS.
Note: If your filter parameters look like https://github.com/mikebrady/shairport-sync/wiki/Digital-Signal-Processing-with-Shairport-Sync I advise some caution: You have a lot of VERY high Q features are 6 Hz/12 Hz (very low freqencies). These could be numerically unstable and you really don't need them. I would eliminate all sections below 30 Hz and just make sure you have decent highpass filter somewhere in the signal chain.


time downsampling vs. frequency downsampling $x[n]_M$ is a finite length sequence of length M.

if:

$$
y = x[nN]_M    \tag{1}
$$

is called downsampling in the time-domain.

then what do you call the process of converting going from a M-point to an N-point DFT?

$$
y[n] = \Bigg[ \sum_{k=-\infty}^{\infty} x[n-Nk]_M \Bigg] R_N[n] \tag{2}
$$

does that have a name?  like frequency downsampling? frequency rebinning...frequency re-pie cutting or something like that?
just guessing...

($R_N[n]$ = window function, which is 1 between 0 and N-1 and 0 otherwise.)

$$
Y[k]_N = X(e^{j\omega})\bigg|_{\omega=2\pi k/N} \tag{3}
$$

rebinning $Y[k]_N$ to $Y[k]_M$ for example.

>

Your second equation suggests that $y[n]$ is obtained by performing an $N$-point inverse DFT of an $M$-point DFT $X[k]$ of an $M$-point sequence $x[n]$.

The result can be interpreted as shifted copies of $x[n]$ are superimposed, causing possible aliasing in time, which depends on whether $ N \geq M$ or $N < M$.

If $N \geq M$, there is no aliasing and $y[n]$ will be an $N$ point sequence whose first $M$ samples will be identical to $x[n]$ for $n=0,1,...,M-1$ and remaining $N-M$ samples will be zero-padded.

$$
y[n] = \begin{cases} { x[n] ~~~,~~~ n = 0,1,...,M-1 \\ ~~~~~0 ~~~,~~~ n= M,M+1,...N-1} \end{cases}
$$

Otherwise if $N < M$, then $y[n]$ will be a time-aliased version of $x[n]$, as shifted copies of $x[n]$ would be overlapping each other. If you want any non-aliased samples in $y[n]$, then $N > M/2$ in general and in such a case $y[n]$ could be defined as:

$$
y[n] = \begin{cases} { x[n]+x[n+N] ~~~, n = 0,1,...,M-N \\ x[n] ~~~~~~~~~~~~~~~~~~~~~,~~~ n= M-N+1,...,N-1 } \end{cases}
$$


Non-Uniform Quantization I was reading a research paper on companding schemes for non uniform quantisation. In paper one of the motivations for non-uniform quantisation is that distortion at larger amplitude values is less audible than distortion in lower amplitude values but I can find why this happens?

The idea behind nonuniform quantization (for example an a, or mu law for telephony speech) is to make better use of the input signal amplitude distribution in lossy source coding applications..
Over a period of time, it's emprically observed that, most of the speech occurs at low amplitude levels, and occasionally at high amplitude (full dynamic range) levels.
When this is the case, a uniform quantizer would be wasting its bits (its quantization levels) into many of those statistically less used louder signal levels.
In order to make efficient use of the quantizer (number of levels and bits), and to minimize the average quantization error (the distortion), it will be wiser to allocate more bits (finer levels) to the (often used) lower signal levels, while assigning less bits, and (coarse levels) to the (less used) louder signal levels.
The technical term for such a modification to a uniform quantizer is called a pdf optimized Max-Lloyd quantizer. The compander is a simplified example by assuming a fixed typical pdf of input speech and assigning a logarithmic non-uniformity to the step-sizes. This way, softer passages will be assigned more bits and finer resolution while louder passages will be assigned less (than uniform case) bits.
Futhermore, the distortion at the louder passages is perceptually less audible, due to psycho-acoustic principles, in particular to the masking effect; presence of a loud tone, suppresses the perception of nearby quiter tones, and makes them inaudible. Hence the quantization noise will be less audible while occuring with louder tones.


Retrieve pixel values from the output jpeg image I have this jpeg image (abc.jpg) and I change some pixel values of this image and store this on hard disk, to later see the changed the pixel values. The name of the output jpeg image is say, xyz.jpg. The problem here is in the image xyz.jpg. The changed pixels does not give me the exact values which I changed them to.

Example:

In abc.jpg at the pixel position (x,y)=(1,10) the pixel value is (231, 120, 120). Now, I want to change it to (230, 120, 120). Now, I save this output image as xyz.jpg. When I try to retrieve the values from position (1,10) I get the pixel value as (232, 120, 120). Now you see the problem, I can't get the value to which I have changed. I am doing all this in MATLAB. And yes, I have tried to save the output image as .tiff, .bmp, .gif. But the size of these formats is very large as compared to .jpg. Please suggest some solution to this problem.

Ask yourself: why do you think the TIFF and BMP files are so large while JPEG files are small?

Answer: JPEG is a lossy compression format and achieves small file size by removing details from your image imperceptible to the human eye. By design, JPEG cannot guarantee that the decompressed image will be pixel-by-pixel identical to the compressed image - it just preserves the overall content the human eye can perceive. You could try increasing the quality parameter, but this probably won't solve the problem in your case. Brightness levels will still be quantized, and color information downsampled.

Some "solutions":

Do not use JPEG, use a lossless compression format.
Do your image editing/modifications directly in the compressed domain. I doubt there are tools in Matlab to do that, and compressed-domain image modification is not as straightforward as writing into an array of pixels!


How should I downsample my sensor? I have 2 sensors, sensor A is sampling a signal at 1kHz, and sensor B is sampling the a signal at 100Hz.

Every timestamp that sensor B takes a measurement at is contained in sensor A. That is, sensor B takes a measurement essentially every 10th measurement taken by sensor A.

How should I downsample sensor A so that both data streams are 100Hz. Can I just pluck every 10th sample from sensor A, or should I low pass filter sensor A first? What frequency should the low pass filter be set at?

EDIT: It seems to me there should be no difference between taking every 10th sample or lowpass filtering if the signal I care about is less than 50Hz. Is that true?

If you can assure that the physical signal being measured by A has no detectable components above 50Hz, then you could take 1 every 10 samples.

However, if you have CPU power available, it is much better to filter the signal (lowpass, with stop-band starting below 50Hz) before downsampling. This way, components above 50Hz won't be aliased down to your working bandwidth. Also, this process will reduce quantization noise from the ADC, which is a good thing.


Java amplitude processing I have some code that collects amplitudes from a complex sound like a song. I now want to plot a graph of amplitude vs time, for the graph plotting I use JavaFX. However I get about 9 million amplitude shorts for each channel, so how to I select a few samples to plot (as plotting 9 million points is obviously out of the question)? Do I need to use the sample rate or sample size? Is there a particular way I should go about this?

Is this post more suited for StackOverflow?

AudioInputStream in = AudioSystem.getAudioInputStream(file);
AudioInputStream din = null;
 AudioFormat baseFormat = in.getFormat();
  AudioFormat decodedFormat = new AudioFormat(AudioFormat.Encoding.PCM_SIGNED,
  baseFormat.getSampleRate(), //sample rate
  16,                         //sample size (bits)
  baseFormat.getChannels(),   //channels
  baseFormat.getChannels() * 2,//frame size
  baseFormat.getSampleRate(), //frame rate
  false);                     //little-endian
  din = AudioSystem.getAudioInputStream(decodedFormat, in);
  byte[] temp = new byte[4];

  //....
  //leftDos and rightDos are streams of shorts
  //....

  int i = 0;
  while (din.read(temp, 0, 4) != -1) {
      if (decodedFormat.getChannels() == 2) {
          leftDos.writeShort(temp[1] * 256 + temp[0]);
          rightDos.writeShort(temp[3] * 256 + temp[2]);
          i++;
      }
  }

 //...retrieve shorts...//
 short[] left = new short[(bytes.length / 2)];

 int len = 0;
 try {
     while (1 > 0) {
         left[len] = leftDIS.readShort();
         len++;
     }
 } catch (EOFException e) {
      System.err.println("End of stream");
 }

There are a number of ways to do this, but I like this way:

// very rough pseudocode
step_size = 9000000 / width_of_window;  // 9000000 or however many points
for each pixel in width of window
  figure out what portion of the audio stream this x value corresponds to
  take step_size shorts corresponding to position in audio stream
  find the minimum and maximum values in that portion of the stream
  draw a vertical line at the current x value with position and length corresponding
      to the min and max values you found

So basically you're displaying the min and max values of successive portions of the stream in each column of your draw box, or canvas, or whatever you're using.


How to compute transfer function $G(s) = \exp \left( - \sqrt{s} \right)$ in Matlab / Simulink? How to compute transfer function
$$G(s) = \exp \left( - \sqrt{s} \right)$$
in Matlab / Simulink?
I am trying to calculate a PID controller for this function. This function describes heat transfer via conduction. I want so see first the step time response and impulse response for this system. I usually use for other functions the tf function from Matlab, including IODELAY, but in this case i can't compute $\sqrt{s})$.The error is:Undefined function 'sqrt' for input arguments of type 'tf'. Is there any other way to compute this transfer function?

Hi again, i just edited it because i found more information.In my case i guest that T=1.

For Figure 2.21 i have this:

The inverse Laplace transform of $G(s)=e^{-\sqrt{s}}$ can be computed analytically:
$$g(t)=\mathcal{L}^{-1}\left\{e^{-\sqrt{s}}\right\}=\frac{1}{2t\sqrt{\pi t}}e^{-\frac{1}{4t}},\qquad t>0\tag{1}$$
[I've cross-checked this result by looking it up in one of my old math books.]
So now you have an analytical expression for the impulse response $g(t)$.
The step response can also be computed analytically. It is given by
$$s(t)=u(t)\int_0^t g(\tau)d\tau=\textrm{erfc}\left(\frac{1}{2\sqrt{t}}\right)u(t)\tag{2}$$
where $u(t)$ is the unit step function, and $\textrm{erfc}(\cdot)$ is the complementary error function.
The result $(2)$ was obtained with a little help from WolframAlpha.
And that's how the impulse response and step response look like:


Convolution theorem for "weighted" convolution? Consider what I call a "weighted" convolution of a two-dimensional signal (image) with itself:

$f({\bf r}) = \int d {\bf r}' \,  g({\bf r}') s({\bf r},{\bf r}') g({\bf r}-{\bf r}')$

where $s$ is a weighting function, with $s({\bf r},{\bf r}') = s({\bf r}+{\bf r}',-{\bf r}')$. My ultimate goal is to perform edge detection on $g$, i.e., find discontinuities in $g$, from knowing only $f$.

If $s=1$, I believe I can obtain $g$ from $f$ using the convolution theorem and a deconvolution scheme, but with a general $0<s<1$ this is no longer possible.

Questions:

Does the integral I write above have a formal name?
Is there something similar to the convolution theorem for such a "weighted" convolution?
Is there any other result that may help me perform edge detection on $g$ from knowing only $f$?

1) The terms autoconvolution or selfconvolution can be found in the literature, without the kernel $s$. So you would probably have a weighed, kernel or windowed self-convolution.
However, i foresee two potential caveats on your formulation. First, i do not understand why you require a two-dimensional periodic signal (and nothing on the kernel $s$?). This could result in infinite values of your integrals. Second, don't you assume some kind of symmetry in $s(r,r′)$?

2) Suppose you can solve your system with $s=1$, then with appropriate symmetry on $s$ you may be able to factorize $s$ and slit it over $g(r')$ and $g(r-r')$. Unfortunately the convolution theorem is likely to give you a family of solutions $g$ knowing $f$ only.

3) You may find useful to explore a multiscale form of your weighted self-correlation operator.

Affine invariant pattern recognition using Multiscale Autoconvolution
Rahtu E., Salo M., Heikkilä J.
IEEE Trans Pattern Anal Mach Intell. 2005 Jun;27(6):908-18.

An iterative deautoconvolution algorithm for nonnegative functions
Kerkil Choi and Aaron D Lanterman
2005 Inverse Problems 21 981

Autoconvolution and panorama: Augmenting second-order signal analysis
Douglas, S.C. and  Mandic, D.P., 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)


How to discard coefficients in the frequency domain without losing too much quality? I was reading a text that shows that one can throw away 66% of the coefficients and gets an image that is still acceptable.

I tried to replicate it using python (last example) and selectively throwing away some of the coefficients but I couldn't get a great image, mine is looking like a ghost.

How can I discard that not so important to the image?

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np
from scipy import signal

# loading the image
img = mpimg.imread('i/super_mario_head.png')

# image in frequency domain
fft_original = np.fft.rfft2(img)

# discarding some coefficients (except initial square 80x80)
for x in range(300):
    for y in range(300):
        discard_x = x > 80
        discard_y = y > 80

        if discard_x and discard_y:
            fft_original[x,y] = 0

# reverse from frequency to spatial
i_fft_original = np.fft.irfft2(fft_original)

f, (plt1, plt2) = plt.subplots(1, 2, figsize=(15,15))

plt1.axis('off');plt1.set_title('Original');plt1.imshow(img[:,:,0], cmap='gray');
plt2.axis('off');plt2.set_title('Reversed');plt2.imshow(i_fft_original[:,:,0], cmap='gray');

You are throwing away coefficients before checking to see if their current magnitudes are small (smaller than typically visible).

(You might also want to rescale what's left to preserve total energy.)


Can wavelets be used for complex texture discrimination? I've recently been studying wavelet analysis with a view to differentiating certain areas of texture images where the texture differs from the background pattern (which is quite random); for example a small ink smudge on a photo of sand or similar.

I've seen a lot of applications of Gabor wavelets used to pick out textural differences in things like fabrics, and other regular type patterns but I've not been able to find any applications of things like Daubechie/Harr wavelets for anything other than image compression.

Can different wavelet families be used to pick out textural differences such as my example or am I wasting my time?

I think you should just dig deeper and get your hands dirty. I have used Mathematica as it is the fastest way I can model something like this

As an example:

(* img *)

dwd = DiscreteWaveletTransform[img, BiorthogonalSplineWavelet[1, 3], 4]

Here's how the decomposition tree looks like

You can then extract every coefficient

coeffs = dwd[{___, 0 | 1 | 2 | 3}, {"Image", "ImageFunction" -> Identity, ImageSize -> 120}]

From left to right {0}, {1}, {2} and {3}.

And you can compare those to the coefficients after applying the 1 + Log[coeff^2] function to them


Object Detection Algorithms I'm a complete beginner to object recognition.

I need to identify and count how many of a certain type of Lego, in a picture full of different pieces of Lego, using MATLAB.

The pieces of Lego I need to identify, are of a certain shape, colour and size.

I need suggestions for algorithms that can perform this task, and links to where I can find more information about them.

From the research I've done, so many are focused around facial recognition, and provide no help for me!

I followed a tutorial that used SURF points to identify points of interest, in the scene pictures and the picture of the Lego picture to find. But it wasn't able to pick up on all the aspects of the Lego piece, and therefore didn't have a chance of finding it in the scene.

If there aren't any algorithms that can perform this, how can I at least separate the pieces based on their colours?

Any suggestions would be really handy! Thanks in advance.

Edit:

Sample picture where I will need to identify certain pieces:

SURF works well for objects with distinctive textures. Unfortunately, all lego blocks have essentially the same texture, so it is not a good cue for telling them apart.  The other difficulty is that your blocks can have different 3D orientations.

The good news is that they are brightly colored. The Color Thresholder app is a good way to segment out objects by color.  Once you have identified individual blocks, you can try to use other cues, such as size and shape to identify different blocks of the same color.  Try using the regionprops function.  For example Eccentricity can help you distinguish between a square block and an elongated block, assuming they are viewed from the top.

Unfortunately, I can only give you a starting point, not the complete solution. You are going to have to do some experimenting here.


How is non-integer upres downres (1080p game on a 1440p screen etc) calculated Please redirect me if this is not the right place to ask.
I was unable to find a straightforward answer to this question. What I am wondering is, if I am playing a 1080p game on a 1440p display, what values will be actually displayed on the screen?
A very simplistic solution would be some form of interpolation, like bicubic. Seems reasonable enough, but is it actually what they do? What if I am playing a 1080p game on a 2160p screen (factor of 2)? Is there any significant difference between integer and non-integer scaling factors?
What if I am forcing the game to render more pixels than the display can show, which some games allow (run the game at 1440p but connect a 1080p display). You can set the render scale to 133%. If you do something like this, the rendered image will have more pixels and more fidelity, but almost none of them will have a 1:1 correspondence to the display. You could still interpolate the 1080p pixel values based on a 1440p rendered image, but intuitively, this seems like a huge disadvantage, but when you do this in practice, the image indeed looks better than when you have 100% render scale where every displayed pixel has a 1:1 correspondence to the rendered image.
If I rendered 2160p and displayed this on a 1080p screen, you could average every 2x2 pixels and display it on a single pixel, this seems simple enough, but is it actually what is done or is it indeed optimal?
Essentially I am asking if there is any genius trick to upresing and downresing images
Another heuristic question: Which would look better? A native 1080p 27" screen displaying a 1080p game, or a 1440p 27" screen displaying a 1080p game? Surely the native 1080p screen must look better? But by how much?

Assuming modern 3D games here.
most games you are aware of are not fixed-resolution – there's very rarely a 1080p game these days! The GPU only renders things onto a buffer the size of your screen's resolution when asked to (in the last step of the processing chain, basically), so games can simply have a fully relative-coordinate-system going on till the end.
So, the answer is: I bet what you think is a 1080p game isn't. Other than that, the physical dimension of the screen make absolutely no difference to the image generation – it's only the last step, where, no matter what the output resolution is, an operation is applied to get pixel values out of a 3D scene definition, consisting of polygons in a 3D coordinate system that has nothing to do with pixels.
You of course need to use appropriately many polygons for high resolutions; the textures you map onto these (and thus, things like static elements, menus etc) need to be right for the rendered resolution etc.
But: there's nothing "special" about the resolution; in any case, there's a projection from one coordinate system to the screen coordinate system (pixels) in a last (or close to last) step.

You could still interpolate the 1080p pixel values based on a 1440p rendered image, but intuitively, this seems like a huge disadvantage, but when you do this in practice, the image indeed looks better than when you have 100% render scale where every displayed pixel has a 1:1 correspondence to the rendered image.

That very much depends on the game you're looking at. In effect, shrinking with an appropriate scaling algorithm is a form of antialiasing setting. That's probably what you perceive as quality advantage: smoother edges.

If I rendered 2160p and displayed this on a 1080p screen, you could average every 2x2 pixels and display it on a single pixel

This should not be done (and rarely if ever is done today); it's a bad way of downscaling things. Signal theory tells us that when filtering with a boxcar filter (in this case, 2 samples wide) in spatial domain, you get ringing, and things look jagged.
GPUs have scalers, but even as simple shader, better scaling algorithms are easy to code.


Is it possible to just remove the noisy part of signal? Is there a way to scan through a signal and remove the parts that don't look "normal". This would mean using more advanced then a simple pass band filter, but something that can look at a period in the signal and remove it if that period didn't fit into a "template".

The short answer is no, unless you know exactly what the noise is, you can not in any practical sense remove the noise. There are things you can do to maximize SNR but no technique exists to remove all noise.


How to create a control loop using a 100Msps digital signal as input I would like to know what type of microcontroller should I buy for this control loop task.
Scenario:
I have a digital signal coming from the following equipment:
$\boxed{\text{Ultrasonic}\\\text{transducer}} \longrightarrow\boxed{\text{ADC}}\longrightarrow\boxed{\text{DSP Board}\\\text{(some high and low pass filter)}}\longrightarrow\boxed{\text{USB 2.0}\\\text{(40MBPS) output}}$
I want to be able to use this digital signal coming trough USB as input in my controller and calculate the TOF (time-of-flight) of the echo and react with my controller every $1\textrm{ ms}$ or less using a PID control law.

How should I deal with this digital signal output in my microcontroller, to be able to detect the echo in generate a reaction in less than $1\textrm{ ms}$?

Any suggestion of microcontroller that I could use for this task?

And how can I deal with such high speed input?

No Microprocessor is able to deal with such extreme rates as 100MS/s. That is clearly the domain of custom ICs and FPGAs, bringing the data to a very capable CPU, not something small and embedded.

Regarding USB:

I don't know your sampling bit width, but 40MB/s is not very much if you need to transport 100MS/s – you get 12bit per sample, in the ideal, zero-overhead case (which doesn't exist).

As a matter of fact, I know a series of devices very well that can deliver 100MS/s and more to a PC – these work through expensive interface such as a proprietary external PCIexpress interface, or better even through 10Gbit Ethernet (which you'd normally only find in data center / ISP equipment).

With a single 100MS/s stream you can easily flatten a single CPU core of a pretty capable x86_64 computer – really, that MUCH MUCH more data than what "microcontrollers" handle.

This might really be a case of insufficiently discussed system design – Nyquist and my understanding of the term "Ultrasound" say that you'd never need more than twice the highest signal frequency as sampling rate. I wouldn't expect ultrasound sonar to have frequencies higher than maybe 4 MHz, but I might be wrong.


Please, help me interpret my Simulink simulation of AM. Why does carrier not appear in the spectogram? Carrier is 2 * pi * 25, modulating signal is 2 * pi, all the amplitudes are 1s and zero phase.

I realized that I could add a sine wave right at the end. I get a carrier in frequency domain, but my signal in time domain looks strange.

I see how you added the carrier by multipying by $1 + 0.5cos(2\pi f t)$ (or you used sine, wouldn't change it), and that approach seems fine to me.

It looks like you do actually see the carrier in your plot! What I see from your plot does appear to be a signal at +/-25 Hz which is what we would expect to see for $cos(2\pi 25 t)$ (or sine if you used that)

However also note that your plot claims the resolution bandwidth is 4.88 Hz. Your modulating signal of 1 Hz would not be visible unless you decrease that resolution bandwidth! Right now all signals (25 Hz carrier and +/-1 Hz modulation sidebands) are within that bandwidth so would appear as one signal which is what you see.

To increase the resolution bandwidth, you need to increase the total simulation time (which would increase the number of points in the FFT if you keep the same sampling rate, or you could decrease the sampling rate with the same number of samples or do both). Resolution bandwidth is approximately the inverse of the simulation time, but that bandwidth would be increased further than $1/T$ if you are using windowing in the FFT. From the plot it does not look like windowing is used, in which case your simulation duration would be predicted to be $1/4.88$ Hz or $204.9$ ms. With your sampling rate of $5$ KHz this would suggest $1024$ samples. (It may be a sliding FFT block size over 1024 samples as I also see T = 10s). Thankfully your plot tells you what the bandwidth is so just increase the total FFT size, decrease your sampling rate or both --to see your $1$ Hz sidebands clearly you should make the resolution bandwitdth on the order of $0.1$ Hz.


Terminologies - lags, order in time series model I am facing some difficulties with the terminologies - lag ($p$) and sequence length (number of data points) (N) used in time series model such as Moving average and Autoregressive model. Considering a Moving Average model described by the equation

$$y(n) = \sum_{i=1}^{p} h(i) x(i-n)$$
 If the number of bits in the input is N=16 defined by $\mathbf{x} = \{1,0,1,1,0,1,0,1,1,0,1,1,0,0,1,1\}$ and the number of lags are p=10
then

Q1: should the model generate a time series of length 'N=16` i.e, would the output of the above model $\mathbf{y} = [y_1,y_2,\ldots,y_N]$ contain 16 elements where $n = 1,2,\ldots,16$?

I am writing out the full expression. I can write the MA process in vector form as $y_n = \mathbf{h}^T \mathbf{x}_n$ where $\mathbf{x}_n = [x_n,x_{n-1},\ldots,x_{n-p+1}]^T$ is the input signal vector and $\mathbf{h} = [h_0,h_1,\ldots,h_{p-1}]^T$ is the MA system coefficients. For example, in programming using Matlab the model would look like for

h = [1,0.2,0.3]; %channel coefficients
N=16;%number of data points
x = rand(1,N); %input data vector
p=2;
y(1) =0.0;
y(2) = 0.0;
for n = 3:N
y(n) = h(1)*x(n) +h(2)*x(n-1)+h(3)*x(n-2);
end

Since the order p=2 would the physical interpretation be that the channel can take in only 2 data values as input? Also, does p=2 mean that the channel has 2 paths through which the data input travels?

Q2: Is there a relationship between the number of input elements and the number of lags?

Please correct me if I have misinterpreted the concepts. Thank you.

Q1: should the model generate a time series of length 'N=16` i.e, would the output of the above model $\mathbf{y} = [y_1,y_2,\ldots,y_N]$ contain 16 elements where $n = 1,2,\ldots,16$?

If one thinks about your question in terms of FIR filtering, then an input signal, $x$ of length $N$ filtered with an FIR filter of length $p$ ($p$ coefficients) will have an output of length $p + N - 1$.

Note: there are ways to define the output so that only $N$ outputs are obtained, but you need to understand that is what you are doing.

  Since the order p=2 would the physical interpretation  be that the channel can take in only 2 data values as  input?

It means the channel does not change the input much at all. It's a two coefficient FIR filter.

  Also does p=2 mean that the channel has 2 paths through which the data input travels?

One way to interpret the FIR (MA) nature of the filter you are using is that it involves $p$ paths with a time delay of $nT$ (where $T$ is the sample period).

  Q2: Is there a relationship with the number of input elements and the number of lags?

In general, no, there will not be a relationship.

  (1) Is this related to "multipath" often used in Rayleigh fading? If there are $p$ paths then would the transmitted signal travel in two directions?

The usual method of modeling multi-path is to say that the receiver receives:
$$
y^c(t) = \sum_{i=1}^{p} h^c_i x^c(t - d_i)
$$
where this is the continuous-time version so the $d_i$ are not integers and $d_1 = 0$.

Moving to the discrete-time version means you need to sample everything using $t = nT$.  So if you stick with $p$ components, you can only really have $d_i = (i-1)T$.

  (2) What is the general way to set $nT$? Using Matlab, I did not do any thing in my code. Can you say what is $nT$ for my example?

$T$ is determined by your modeling. How did you relate the actual continuous-time problem to the discrete case?


Downsampling audio for use in Machine Learning I'm trying to use the work (Neural Networks) done in this repo:
https://github.com/jtkim-kaist/VAD

It says this:

  Note: To apply this toolkit to other speech data, the speech data
  should be sampled with 16kHz sampling frequency.

I've got speech data at 48khz. I've read in places that reducing sampling rate is a complicated process, you can't just remove every nth datapoint, you have to filter things...

Is this necessary if I only intend to use the data in the Neural Network toolkit provided by the repo I linked?
If so, is there an industry standard method for changing sample rate?

I realise that it probably depends on what features are being used. However the feature that is used is this:

  MRCG (multi resolution cochleagra) concatenates the cochleagram features at multiple spectrotemporal resolutions

This is a ruddy complicated feature! Lets pretend we're just using a Melspectogram (unless you're willing to answer the question from the perspective of MRCG's).

Neural networks are likely to use features of a Melspectogram that we wouldn't think of. This makes me think it is unwise to train the Neural Net using downsampled speech data unless we intend to predict using 48khz data downsampled to 16khz forever after...

What do you think? Can I use my 48khz data - downsampled with no filtering - with the expectation that the model will work for prediction on real 16khz data?

And then for future readers sake, how about the other way? Say I had an 8khz file, could I increase the sample rate to 16khz without filtering?

Is this necessary if I only intend to use the data in the Neural Network toolkit provided by the repo I linked?

Yes.

Whether or not you are downsampling (instead of just decimating) has nothing to do with classification performance but rather, it is to preserve (as much as possible) the information  contained in the signal.

When changing the sample rate of a signal, you need to account for aliasing. Otherwise, the signal will be distorted in such a way that it would be impossible to reconstruct the original signal from the samples.

  If so, is there an industry standard method for changing sample rate?

Yes.

The "right" way to do it when downsampling is to first apply an anti-aliasing filter and then decimate the signal but when upsampling, you first upsample and then apply interpolation (which can also be expressed as a filter). Various platforms provide functions to do just that (e.g. Python, MATLAB).

  I realise that it probably depends on what features are being used.

To some extent. Perhaps you could simply leave everything unchanged, provided that you adjust certain points in the algorithm that are specifically hardwired for a sampling frequency that is different than the one you are using. However, the bandwidth for speech is limited. If your classifier is also trained on what is NOT speech, having more bandwidth means that you can discriminate more things. If the majority of your dataset is pure speech then having more bandwidth would not make a huge difference.

  Can I use my 48khz data - downsampled with no filtering - with the expectation that the model will work for prediction on real 16khz data?

No.

If you downsample without filtering you are very likely to produce distorted signals.

  And then for future readers sake, how about the other way? Say I had an 8khz file, could I increase the sample rate to 16khz without filtering?

No.

When you increase the sample rate, you have a more serious problem. You now have to fill "gaps" for which you do not have data. The filtering (or interpolation) stage fills in those gaps but does not add any new information. When you downsample, you are throwing away information. Upsampling cannot replace it.

Perhaps it would be beneficial for you to look up some basics of DSP such as sampling, sampling and reconstruction, aliasing and anti-aliasing filters. Also, to further understand what your features look like, how do they work and how to choose and use them, you might want to look up more information on the Short Time Fourier Transform and Multirate Signal Processing (more generally).

Hope this helps.


Behaviour of microphone with high volume input I have a microphone which produces unexpected results when exposed to very loud input. The mic is a cheap electret headset mic and I am recording into a DAW. For very loud, noisy signals (eg. blowing on the mic) I get an unexpected waveform. In the MATLAB plot below, you can see a snippet of about 0.5 seconds of audio which demonstrate what I'm talking about. The signal is clearly clipping at -1.0, but is not getting near the other end of its range at 1.0.

In the same signal, the clipping behaviour can then reverse and start clipping almost constantly at 1.0 and not reaching -1.0. It's as though there has been some sort of wandering, fluctuating DC offset introduced. Can anyone shed some light on what is happening?

An overview of the analog signal path in a system such as yours:

Clipping can take place:

In the microphone as mechanical soft clipping, distorting the waveform rather than clipping exactly at some value.
In the built-in impedance converter field effect transistor (FET) of the electret microphone capsule. (I think this is very much like hard clipping)
In a separate preamplifier as analog hard clipping.
At the analog-to-digital converter of the audio interface, as hard digital clipping.
On the digital signal if any digital processing is done on integer data, as hard digital clipping.

Because the clipping is hard at numerical values -1 and +1, digital clipping must be taking place. Because it looks like there is hard clipping also at other numerical values than -1 and +1, there is probably hard clipping taking place at more than one stages separated by an alternating current (AC) coupling. Blowing directly into the microphone is quite an extreme disturbance, so there is probably also the 1st and the 2nd kind of clipping going on.

There are a few places in your system where there can be an AC coupling or equivalent:

In the microphone, as the pressure in the back of the omnidirectional closed-back capsule slowly equalizes with that of the outside world due to a (by design) leaky seal.
On the analog signal. AC coupling may be built in into the plug-in power supply to the microphone, into the preamp and into the audio interface.
Digitally in the audio interface or in the digital audio workstation (DAW). Check the settings!

If you blow into the microphone long enough then the AC coupling begins to see the average pressure level of that as direct current (DC), centering the signal to zero during the disturbance. When you stop blowing it takes some time for the coupling to adapt to the normal atmospheric pressure.

If you want the electronics to recover faster after a disturbance, try adding a capacitor in series with the signal output of the plug-in power supply as yet another stage of AC coupling. The capacitor should be small enough to make the time constant of the resulting resistor–capacitor (RC) filter short enough.


Variance of filtered white noise I was asked a question, as posted here, and the answer given is (A) i.e. $\frac{3}{2} A^2 N_0$.
My solution steps was:

Finding the mean of the output process: Since input is gaussian the output will be also Gaussian, as per property of Gaussian variable.

Finding the mean square value (MSV) of the process i.e. $\mathbb{E}(Y^2)$.

Variance = MSV$-\mathbb{E}(Y)^2$

But I am confused do mean of output process will be zero? Since the input is zero-mean white Gaussian noise if yes how?
Is there any other way to solve the question?

The output process is clearly zero-mean because the LTI system cannot add a mean to the zero-mean input process. The variance of the filtered process is given by
$$E\{Y^2\}=\frac{N_0}{2}\int_{-\infty}^{\infty}|h(t)|^2dt\tag{1}$$
which results in $3A^2N_0/2$.


Fourier Transform of Alternating Periodic Rectangular Pulse I'm having trouble determining Fourier transform of signal. I have 2 ideas on how to solve this problem.
Given the signal is periodic I could use formula for Fourier transform of periodic signals:

$$X(j\omega) = \sum\limits_{k=-\infty}^{\infty} C_k \cdot 2\pi \delta(\omega-k\omega_{0})$$
where $\omega_0 = \frac{2 \pi}{4T}$.

Also, I could make signal $g(t)$ non periodic rectangular signal, and $h(t)$ alternating pulse train, and then $g(t)*h(t)$ should give me signal I needed.
Although I have some ideas I'm stuck solving this problem. Please help. Thanks

The answer is simple.
I will give 3 points to solve it:

The Fourier transform is linear. Hence $ \mathcal{ F } \left\{ \alpha f \left( x \right) + \beta g \left( x \right) \right\} = \alpha \mathcal{ F } \left\{ f \left( x \right) \right\} + \beta \mathcal{ F } \left\{ g \left( x \right) \right\} $.
Shift in time $ f \left( x - {x}_{0} \right) $ equals multiplication by $ {e}^{-j \omega {x}_{0}} $ in Fourier domain.
Instead of solving for the case above, think of the case you have 2 rectangular signals with twice the period with one multiplied by $ -1 $ and shifted.

In your case, just have a look on the signal of the Positive Pulses. It has a period of $ 4T $.
You have another signal. The signal of the negative pulses.
It is basically the same signal as the positive one (It also has a period of $ 4T $) yet it is shifted by $ 2T $ and as it is multiplied by $ -1 $.
So if you know the transform of the positive one, follow my above points and you have the transform of the negative one and their sum.


Difference between 'spatially' and 'instantaneously' independent For a multivariate time series

$$\underline{y}(t)=[y_1(t) \dots y_n(t)]^T$$

What is the difference between the components $y_i(t)$ being spatially independent and instantaneously independent?

The components are instantaneously uncorrelated if
$$
E[y_i(t_1) y_j(t_2)] = 0
$$
for all $i \not= j$ and for all $t_1,t_2$ and when $i = j$ :
$$
E[y_i(t_1) y_i(t_2)] = 0
$$
for all $t_1 \not= t_2$.

The components are spatially uncorrelated if
$$
E[y_i(t) y_j(t)] = 0
$$
for all $i \not= j$ (assuming a change in index corresponds to a spatial change).

Note that correlation is not quite the same thing as statistical independence.


Causal and non-causal systems in real world Causal systems are defined as the systems whose :

Output depends on present as well as past inputs and the impulse response; and sometimes depends also on the past outputs (in case of recursive systems).
Impulse response $h(n)$ consists of only positive samples.

Also, I am trying to link this concept to real world scenario using FIR filter whose design equation is:$$y(n)=h(0)x(n)+h(1)x(n-1)+\ldots+h(M)x(n-M)$$

Based on all this theory, I have questions as given below :

a. I wish to ask if it is possible for any system (non-causal/anti-causal real world system) to have negative samples?? Is it that be it causal or non-causal, all real world systems always have positive samples only? Negative samples are of consideration only in theory?
b. We show past inputs using delay blocks $\left(z^{-1}\right)$. So, this means delayed signals give us past inputs in real world? I am confused. I thought the filter uses past inputs saved in its memory. How can we get past inputs by adding a delay?
c. If I take any real world filter, and suppose I give an input sine wave of some frequency to it. Then, filter decides if this frequency should be allowed to pass or not depending on the past input samples of this sine wave which it had stored in its memory? So, this way the output depends on past inputs as well and so we say it is causal system?

In discrete-time systems, causality is a requirement only when processing (filtering) signals in real time; i.e. when the index $n$ corresponds to a physical time $n \times T_s$. In this case, a non-causal system is impossible to implement since calculating the current output would require inputs from a future time.

In practice, however, many systems are not real time.

Image processing, for example, has indexes for horizontal and vertical coordinates. So it's OK to apply a non causal filter to an image (usual examples being 0-centered FIR filters for blurring or sharpening, or IIR filters applied in both directions to obtain 0 phase distorsion).

In time sequences, many times the signal is stored and then processed in non real time, so you may apply non-causal filters there also.


Calculating output of a system given impulse response and input I'm trying to calculate the output given an impulse response $h(t)=e^{-t}u(t)$ and input $x(t)=\cos(2\pi t)$.

I know I need to use the convolution but given that cosine is periodic, I don't see how I would get an answer to converge. Am I thinking about this the right way or missing something?

Eventhough you could solve this problem using other means, such as frequency domain methods, you could also follow a direct time domain path as the following.
Given the impulse response $h(t) = e^{-t} u(t)$ of an LTI system and the applied excitation input $x(t) = \cos(2\pi t)$ (which is of infinite extend from $t=-\infty$ to $t=\infty$), the output can be written as the convolution integral:
$$y(t) = x(t)\star h(t) = \int_{-\infty}^{\infty} h(t-\tau) x(\tau)d\tau = \int_{-\infty}^{\infty} e^{-(t-\tau)}u(t-\tau) \cos(2\pi \tau)d\tau$$
$$y(t)  = e^{-t} \int_{-\infty}^{t} e^{\tau} \cos(2\pi \tau)d\tau$$
At this point we have to evaluate the integral which can be found on most calculus texts or in an integral table. I will use MATLAB to evaluate it which results in
$$\boxed{ y(t) =  \frac{\cos(2\pi t) }{1 + 4\pi^2} + \frac{2\pi \sin(2\pi t) }{1 + 4\pi^2} }$$
You can simplify the result using $a\cos(wt) + b\sin(wt) = \sqrt{a^2+b^2} \cos(wt - \tan^{-1}(b/a) )$ into
$$\boxed{ y(t) =  \frac{1 }{\sqrt{1 + 4\pi^2}} \cos(2\pi t)}$$
as the solution. The solution has no transient part as the input was applied beginning from $t=-\infty$.


Variance of periodogram estimate of the power spectrum I have been reading chapter 13.4. ("Power Spectrum Estimation Using the FFT") of the Numerical Recipies Book.

Some things related to the expectation value of the "periodogram estimate of the power spectrum" became not clear to me though.

Background

Suppose we have an equally-spaced $N$-point sample  $ c_0 ... c_{N-1} $ of the  discrete signal $ c(t)  $ and its Fourier transform $$ C_k = \sum_{j=0}^{N−1} =
c_j \exp{
2πijk/N}
 $$ $$k = 0, . . . , N − 1 $$
Then we can define the "periodogram estimate of the power spectrum" $P(f_k)$ as
$$ P(f_k) = \frac{1}{N^2}
[{|C_k|^2 + |C_{N−k}|^2]}
$$

with $k = 1, 2, . . . , \frac{N}{2}− 1$ and $f_k = \frac{k}{N\Delta}$ for $k = 0, 1, . . . ,\frac{N}{2}$ ($\Delta$ is the sampling interval).

So far everything is clear.

Now the book says

  "... the variance of the periodogram estimate at a frequency $f_k$ is
  always equal to the square of its expectation value at that frequency.
  In other words, the standard deviation is always 100 percent of the
  value, independent of N!"

My questions

What is the expectation value of the periodogram estimate at a frequency $f_k$? Why (can you show me how to derive it)?
The expectation value doesn't change with $N$, so changing the length of the input doesn't improve the power spectrum estimation. But if I have a non-integral number of periods in my sample, more samples (longer sampling time) means less weight on the end of the signal (where there is a disconnect) so less smearing of frequency components measured due to the disconnect between the beginning and the end of the signal. So in this case a higher $N$ means a better estimation? Where am I wrong?
If the SD is always 100 percent of the expectation value, doesn't that make this simple periodogram estimation (using just a single DFT of $N$ samples) pretty much useless? Should one always opt for Welch's or Bartlett's method e.g.?

Thank you very much in advance!

There are two sources of error in estimating a quantity: variance and bias.  In your question you have only considered variance.

The improvement in the spectrum estimate using the periodogram comes because there is less smearing of frequency components as you include more samples in your estimate of the power spectrum.  Using a longer data window causes less power from one frequency to be leaked into other frequencies.  However, this only reduces the bias of the estimate.

So, as you use more samples the bias is decreased, which means you do have a better estimate.  However, the variance remains unchanged, as stated in Numerical Recipes.  A discussion of this is given on this web page:

http://matlab.izmiran.ru/help/toolbox/signal/spectra7.html


How to plot the Hilbert Spectrum in Hilbert-Huang transform? I'm working on integral transforms for a science initiation project and i got to do a work on the Hilbert-Huang transform.

I've seen on the internet that, after one found the IMFs with the EMD method, you can do a Hilbert Spectrum plot of the absolute value of the componentes $H(\omega,t)$ where $t$ is for time, and $\omega(t)$ is the instantaneous frequency.

Suppose i got 5 IMFs stored in 5 vectors in GNU/Octave (or Matlab, wharever). How can i do such plot?

Final step is pretty straightforward. All you need to do is to apply the Hilbert Transform to each IMF and extract the instantaneous frequency from analytical signal. Instantaneous frequency is given by:
$$\omega(t)=\dfrac{d\phi(t)}{dt} $$
where $\phi(t)=\mathrm{arg}[x_a(t)]$ (unwrapped phase of the analytical signal). Keep in mind that MATLAB (Octave) hilbert function already returns analytical signal.

In fact you might want to re-use this code from Mathworks File Exchange: plot_hht it is a good starting point.


Acoustic measure of a specific room I would like to make an acoustic measure of one room, so I can compare it with a acoustic measure of some other room. Then compare these data. So when I take an acoustic measure of the first room again, I could recognize that i am in the first room again. Next downside is that I would try to use my iPhone's microphone and the build speakers for that task.
Is this even remotely possible?

Can please someone point me in the right direction.

Tricky problem. Any speaker to mic path in a room is actually a linear time invariant system so they can be fully described by impulse response or transfer function. However the transfer functions are very complicated (1000s of degrees of freedom) and vary tremendously with position in the room. The room is also a very effective de-correlator. Signals from the same speaker recorded at two different location will seem mostly uncorrelated.

You could try to use some statistical room acoustical properties: diffuse field power, reverb time, early decay time, reflection density, mean free path lengths etc. However, these are difficult to determine and not that much different between most rooms in the same house.

So unless you can deploy some significant constraints (control the excitation, control location of source and microphone), this will be very hard to do.


Find common features across images Consider multiple images that share a common feature. More specifically, they are the sum of multiple features, and I know that one of those is present in all. What is the best method to extract this common ground?
My naive approach would be to take the minimum over all images for every pixel, but I wonder whether it can be done better.

I am mostly working in Fourier space, but this being a local problem, I expect it to be solved better in real space? A solution in Fourier space would be ideal, though.

For an illustration, see the following images:

What if Goofy was translated and/or rotated?

If by feature you mean a group of pixels (like Goofy) I suggest to try SIFT + SVD.
 (http://weitz.de/sift/)


Harmonic Product Spectrum from FFT output to the fundamental frequency I'm trying to do the HPS algorithm on my FFT output, but I'm stuck here.
My signal is from a bass (cable), 16bits, 48.000Hz.

I split in time domain the array of the whole record into multiple arrays of 24.000 frames (miniBuffer[24000]).

I'm applying a Hamming window once it's split, like this:

miniBuffer [m] = miniBuffer [m] * 0.5 * (1.0 - Math.cos(2.0*Math.PI*m/24000));

Then, I compute FFT on each array.

Finally I try the HPS:

 double mag[] = new double[miniBuffer.length]; //Saves magnitude of the fft values

 double m_max0 = miniBuffer[0];

                    for (int i=0; i < miniBuffer.length / 2 - 1; i++) {
                        double real = miniBuffer[2*i];
                        double imag = miniBuffer[2*i+1];
                        //Gets magnitude
                        mag[i] = (float)Math.sqrt(real*real + imag * imag);

                        if (miniBuffer[i] > m_max0 && i>0) {
                         m_max0 = miniBuffer[i];
                        }
                    }

                 int i_max = 0;
                 double m_max = miniBuffer[0]; //max
                 m_max = m_max0;

                 double sum [] = new double[miniBuffer.length];
                 int fund_freq = 0;

                 //HPS

                 for(int k = 0; k < 24000/48 ; k++)
                 {
                     sum[k] = miniBuffer[k] * miniBuffer[2*k] * miniBuffer[3*k];
                     // find fundamental frequency (maximum value in plot)
                     if( sum[k] > m_max && k > 0 )
                     {
                         m_max = sum[k];
                         i_max = k;
                     }
                  }
                  fund_freq = i_max * 48000 / 24000;
                  System.out.println(fund_freq);

The values of the max magnitude are the frequency that is recording the program, but when i apply the HPS, for example, to a 220Hz sound, the HPS returns me 26 or 30.

I think that this values come from the frequencies limiter that doesn't works (I want to cap before 37Hz and after 440):

                 double hz = ((double)( miniBuffer.length * 2)) /48000;

                 for(int z = 0; z<(37*hz); z++) {
                     miniBuffer[z] = Double.NEGATIVE_INFINITY;
                 }
                 for(int z = 0; z>(440*hz); z++) {
                     miniBuffer[z] = Double.NEGATIVE_INFINITY;
                 }

If I delete this, HPS returns me 0 in every output.

So:

Is the window correctly applied?
Is the HPS a good algorithm for bass pitch detection? (Low Frequency)
What could i use for limit my values in a right way?
How to fix my HPS algorithm?
Some other advice?

Thank you!

If you zero out most or all of the harmonics (you set them all to -inf in your cap), then the HPS algorithm can’t find them.


Vehicle segmentation and tracking I've been working on a project for some time, to detect and track vehicles in video captured from UAV's, currently I am using an SVM trained on bag-of-feature representations of local features extracted from vehicle and background images. I am then using a sliding window detection approach to try and localise vehicles in the images, which I would then like to track. The problem is that this approach is far to slow and my detector isn't as reliable as I would like so Im getting quite a few false positives.

So I have been considering attempting to segment the cars from the background to find the approximate position so to reduce the search space before applying my classifier, but I am not sure how to go about this, and was hoping someone could help?

Additionally, I have been reading about motion segmentation with layers, using optical flow to segment the frame by flow model, does anyone have any experience with this method, if so could you offer some input to as whether you think this method would be applicable for my problem.

UPDATE: I posted this question on stack overflow as well, and had an excellent answer, I've implemented this idea already and it is working surprising well and im now investigating using optical flow in addition to this technique.

Below is two frames from a sample video

frame 0:

frame 5:

Alas, optical flow is a difficult problem too ;-)

Well, to be more constructive, here are a few algorithms that should be worth trying (or have been tried on this particular sequence) :

re-train your bags of features on a databse of vehicles more representative (in size and orientation) to your actual problem in order to obtain better results
use the fact that the ground is a flat plane to do some parametric optical flow (search for an affine flow) or to compute some affine registration between the frames of the sequence. The moving vehicles will then be outlier from this dominant motion
use some optical flow algorithm to compute the flow, then try to classifiy / cluster the optical flow vectors (this is still a widely open problem!). Depending on the language you use, you can use OpenCV's optical flow, the one from TU Graz, D. Sun's optical flow, or even mine ;-). Note however that segmenting the flow will be a non-trivial task that you should probably do in two steps: global (dominant) motion etsimation, then small motion detection.


Amplitude and phase spectrum in MATLAB I have a little problem with the subject Theory of signal. We were told to make amplitude and phase spectrum of cosinus in MATLAB, but we don't have results nor materials from our teachers, so we have to study from many different materials and you can imagine what this is like. So I'm a little confused.

However, I did some graphs, that should be what I want, according to instructions, but I have never seen something like that, so I'm not sure, if that's right.

Do these graphs make any sense to you? I used fast Fourier transform (fft(y) in MATLAB) on the signal and then plot the result of transform versus sample frequency and I thought it's possible to get phase of the signal by angle(fft(y)), so I plotted this on the last graph. I'm very unsure especially by the third graph.

I would like to apologize for incorrect notions I might use, 'cause I'm not a native english speaker.

Thanks for any help.

This is the most simple FT that you should understand. The basic rule is: a sine of frequency f0 in time domain is a delta at f0 and -f0 in frequency. You can do this theoretically, as shown in http://mathworld.wolfram.com/FourierTransformCosine.html.

Those peaks in your graph are consistent with this. However, it isn't true that there is a peak in at 90Hz. The peaks should be at 10 Hz and -10 Hz. There are two possible approaches here:

a) Use fftshift to correct this issue. This will take each frequency to its real place.

Fs=100; %sampling frequency
T=1; %signal length
N=T*Fs;     %number of samples
f=-Fs/2:Fs/N:Fs/2-Fs/N; %frequency vector
x=cos(2*pi*10*t);
Xf=fft(x);
Xf=fftshift(Xf);
plot(f,abs(Xf)) %magnitude

b) Assume the signal is real, so the spectrum is symmetric. It is true in this case, but it doesn't have to. You could just plot the first half of your fft, i.e. from 0 Hz to 50 Hz.

By the way, if you plot plot(f,abs(Xf)/N), you will find that the peaks have an amplitude of 0.5A (being A the amplitude of the sinusoid). The phase plot seems OK to me, especially after fftshift-ing it.


Determining the uncertainty of an autocorrelation My problem should probably be built up from the beginning, so lets start there. I performed a certain experiment 25 times. Every time, the experiment consists of 5000 measurements, and each measurement returns either 1 ('yes') or 0 ('no'). I know that each time I measure 1, the probability that this measurement is correct is $F_1$, while if I measure 0 I know that it's correct with probability $F_0$.

Now, subsequently I want to compute the autocorrelation function of these measurements I use the following formula

with a slightly different normalization, but in principle that should not matter, because what I want to do is fit this autocorrelation to an exponential decay and find the decay time. To do so, I thought I had two options: Either I fit the 25 datasets separately, and then find the mean of the decay time $t_1$ and maybe say something about its variance, or I can find the mean autocorrelation of the 25 datasets, and fit that one. I've decided that that is probably the best option, as the individual datasets can be quite noisy with strange fluctuations, while the mean looks much cleaner.

But now I'm stuck wondering what I should be doing to find the uncertainty in the fit. Should I have introduced some sort of uncertainty due to the imperfect measurements ($F_0$ and $F_1$), or should I have introduced some sort of standard error of the mean for the 'mean autocorrelation' data? That last one seems plausible, as of course there is also a standard deviation in the 25 datasets, but then I get a little confused as to how I should do that. The formula for the standard error of the mean seems to be $\frac{\sigma}{\sqrt{n}}$, but this is a bit unfair as in my 5000 measurements, the autocorrelation drops to 0 after around 50 measurements. For the fit I therefore also only use the first ~100 points of the autocorrelation, as there's no real need to fit the next 4000 points that are all pretty much equal to 0. So should I instead only calculate the mean and standard error of the first 100 points? Would that be a fair way of finding the uncertainty? The subsequent fitting method will already give error bars for the $t_1$, but this of course heavily depends on the uncertainty in the data.

A final sort of sidetrail, what kind of quantitative method of finding how good I fit should I use for exponential decay? Reduced Chi square, Adjusted R squared, or something else entirely? I suppose this is not intended for dsp though, more for a statistics forum, so feel free to ignore it.

As far as I can see the main question is: How to account for the uncertainty in the measurements of $X_t$ (i.e. $F$).

One approach to do this would be to run some sort of Monte Carlo simulation on the data. I.e. for each experiment you can resample a large number of times (say 1000 or 10000) where each point has $1-F$ probability of being flipped. By calculating the auto-correlation of each resample you get a distribution of the correlation for each $k$. Then take the 95% (or another value) confidence interval to get an uncertainty in $R(k)$.

It may also be possible to do this analytically but I'm not sure for binary data (my instinct says not or at least not simply).

This uncertainty in $R$ should be accounted for in your fitting of the data when calculating the decay time. This can definitely be done analytically and may be implemented in whatever you are using to do the fit (I think matlab does). Alternatively a similar Monte Carlo simulation could be done but is probably more effort.

If you think all the data follows the same trend (not unreasonable if it is a repeated experiment) you could calculate a mean for each auto-correlation distance as you suggest. This is effectively taking a manual Monte Carlo with 25 resamples. Then you can calculate the standard error as you suggest (note: $n=25$, not $5000$). I would shy away from this approach as it ignores the uncertainty in you initial measurement which you presumably have a good reason for giving.


Is there a method to calculate the maximum gain of a biquad? Given an IIR expressed as a biquad:

$$
H(z) = {b_0 + b_1 z^{-1} + b_2 z^{-2} \over 1 + a_1 z^{-1} + a_2 z^{-2}}
$$

Is there a way to calculate the maximum gain the biquad attains over the range $0$ and $\pi$?

One way is to calculate the turning points of $H(z) \bar H(z)$ and plug those back into $H(z) \bar H(z)$ but that would result in a horrible equation to solve in large powers of $z$ (probably up to order $7$) or differentiation by parts of a quotient based on terms of $\cos$ (depending on the two ways to tackle $H(z) \bar H(z) $).

Has this problem been simplified down to a ready made formula? Or is there an easier way?

in this previous answer, the squared-magnitude of a biquad frequency response is shown to be:

$$\begin{align} f(\phi) \triangleq |H(e^{j\omega})|^2
& = \left| \frac{b_0 + b_1 e^{-j\omega} + b_2 e^{-j2\omega}}{a_0 + a_1 e^{-j\omega} + a_2 e^{-j2\omega}} \right|^2 \\
& = \frac{B_0 + B_1 \phi + B_2 \phi^2}{A_0 + A_1 \phi + A_2 \phi^2} \\ \end{align}$$

where $ \phi \triangleq \sin^2\left(\frac{\omega}{2} \right) $

and

$\begin{align}
B_0 & = \left(\frac{b_0+b_1+b_2}{2}\right)^2 \\
B_1 & = - (4b_0b_2 + b_1(b_0+b_2)) \\
B_2 & = 4b_0b_2 \\
A_0 & = \left(\frac{a_0+a_1+a_2}{2}\right)^2 \\
A_1 & = - (4a_0a_2 + a_1(a_0+a_2)) \\
A_2 & = 4a_0a_2 \\
\end{align}$

maxima and minima occur when the derivative is zero.

the derivative of the above is zero when

$$ \phi = \frac{A_2B_0 - A_0B_2 \pm \sqrt{A_0^2B_2^2 + A_2^2B_0^2 - 2A_0A_2B_0B_2 - A_0A_1B_1B_2 \\ + A_0A_2B_1^2 + A_1^2B_0B_2 - A_1A_2B_0B_1}}{A_1B_2 - A_2B_1} $$

then evaluating $f(\phi)$ at those two values of $\phi$ is another laborious pain in ass.  and it won't be my ass.


digital IQ modulators when it comes to digital I/Q signal power measurement we know that this power is viewed as $I^2+Q^2$ (as its polar form). This is the power as referred to final analogue signal. Yet the IQ modulator does not convert I/Q to polar form, instead only the final I is output to real world say in the RF modulator.
When I model two I/Q mixers it does not show that final power is $I^2 + Q^2$ but rather just $I^2$:

x = randn(1,1024)+j*randn(1,1024); %random I/Q input
sincos1 = exp(j*2*pi*(0:1023)*.01);
y1 = x.*sincos1; % first digital IQ modulator

sincos2 = exp(j*2*pi*(0:1023)*.1);
y2 = real(y1.*sincos2); %rf IQ mixer giving final analogue output

if you check power of y2 you will see it is close to real(y1) rather than complex y1.
I ask this because for all power dependent applications e.g. DPD we assume power is $I^2+Q^2$
what am I missing here?

Thanks

The power of the final real-valued signal $y_2$ is half the power of $y_1$ because you take the real part (and thus reduce power by a factor 2). The real part of $y_1$ is already the sum of the I-component and the Q-component, not only the I-component:

$$Re\{y_1\}=Re\{x\}\cos\omega_1 t - Im\{x\}\sin\omega_1 t= I-Q$$

So the original formula $I^2+Q^2$ is correct, you just need to apply a factor of $1/2$ if you don't compensate for taking the real-part. This is why in some textbooks the real part operation in a bandpass modulator is 'implemented' as $\sqrt{2}Re\{\cdot\}$; this will keep the power unchanged.