texts: Transfer function estimation of a noisy system

texts:
Overall description
I am trying to estimate a filtering system’s transfer function, given its input and output.
This system takes $x$ as input . This signal is low pass filtered and added to a WGN by the system, which gives the output signal $y$.
\begin{align}
x &= \text{BPSK} + b_1 \\
y &= h*x + b_2 \text,
\end{align}
$b_1$ and $b_2$ being white gaussian Noise
Hypothesis

$b_1\ll b_2 \ll x$
$b_1, b_2$ independent

The purpose of my manipulation is to estimate as best as possible the filter’s frequency response amplitude ($|H|$) and compare it with the applied filtering (which is known – red plot on next figures).
I had some trouble with this task and it would be very helpful if someone could help with an additional signal conditioning or identify an algorithmic issue.
Algorithm
Theoretically, a simple estimator based on PSD (Power Spectral Density) of input and output signals can be performed with following equations:
$$H(\omega) = \frac{S_{xy}(\omega)}{S_x(\omega)}$$
with $S_{xy}(\omega) = \text{FFT}\left(R_{xy}(\tau)\right) = \text{FFT}(\mathtt{xcorr}(x,y))$ and $S_{x}(\omega) = \text{FFT}\left(R_{x}(\tau)\right) = \text{FFT}(\mathtt{xcorr}(x))$.
For a better signal to noise ratio, the signals ($x$ and $y$)  are split into bursts to compute a mean FFT.
Results
The output of this estimator is also given hereafter (Estimated TF = estimated transfer function). It can be seen that the estimator is noisy.

In this step I have tried to enhance the estimator by adding median filter on the transfer function estimate. However, this did not improve the estimator.

It seems that the Frequency response amplitude estimate shows  spectral hollows compared to the input filtering, which are inducing an important error on the estimate.
Questions

What kind of filter (or another way to estimate the transfer function) can I use?
Could a signal conditioning on $x$ or $y$ enhance the
estimator accuracy?

The following Matlab function is used to compute the TF estimation. Please note that nb_means parameter is arbitrary equal to 16 in my computation, this parameter can evolve.
function [ H1, G1 ] = estimateTF_DSP( sigX, sigY, nb_means )
if nb_means > 1
    sigX_I_reshaped = reshape(real(sigX),length(real(sigX))/nb_means,nb_means)';
    sigX_Q_reshaped = reshape(imag(sigX),length(imag(sigX))/nb_means,nb_means)';
    sigY_I_reshaped = reshape(real(sigY),length(real(sigY))/nb_means,nb_means)';
    sigY_Q_reshaped = reshape(imag(sigY),length(imag(sigY))/nb_means,nb_means)';

    for k=1:size(sigX_I_reshaped,1);
        sig_X = sigX_I_reshaped(k,:)+1i*sigX_Q_reshaped(k,:);
        sig_Y = sigY_I_reshaped(k,:)+1i*sigY_Q_reshaped(k,:);
        % XY
        Rxy = xcorr(sig_X,sig_Y);
        Sxy(k,:) = fft(Rxy);
        % X
        Rx = xcorr(sig_X);
        Sx(k,:) = fft(Rx);
    end

    Sx_m = mean(Sx);
    Sxy_m = mean(Sxy);

else
    % XY
    Rxy = xcorr(sigX,sigY);
    Sxy_m = fft(Rxy);
    % X
    Rx = xcorr(sigX);
    Sx_m = fft(Rx);
end

H1 = (Sxy_m./Sx_m)./(2*pi);
G1 = 20*log10(abs(H1));
end


texts: The DC bin of an FFT should be the mean of all the other bins in the spectrum, why does the FFT produce 0?

texts:
I am experimenting a bit with FFT's. From what I have read the first bin of an FFT is called the DC bin -- that is it is the mean of the other components in the FFT. I have also read that it should be close to 0.

When I actually perform the FFT and take the absolute value, what I find is that the first element of this is always 0. This makes sense mathematically if you consider the definition of the DFT itself.

Is there a best practice to reconstructing this first DC bin?

So it seems that because I was making my test FFT with a sine wave as the transient signal, this is what caused the 0 valued first FFT bin. I changed it to a cosine and got a non-zero first element...And the same with any random data. I should have thought a bit first.


texts: Why does Simulink generate this code for a PID controller?

texts:
For the Simulink PID Controller model

The Simulink generated code (rewrite for better understanding) is:
#define PERIOD 0.005
double PID_step(double err,double P,double I,double D,double N){
  static double Integrator,Filter;
  double POut,IOut,DOut,Out;

  POut= err * P;
  DOut = ((err * D) - Filter) * N;//LineA
  Out = (POut + Integrator) + DOut;
  Integrator += (err * I) * PERIOD;
  Filter += PERIOD * DOut;//LineB
  return Out;
}

Line A and B are the code for the differential part, err is the input and DOut is the output. The Z-domain transfer function is $$\frac{DN}{1+\frac{NT_s}{z-1}}.$$
I want to derive the code from the transfer function.
Let $$\frac{Y(z)}{X(z)}=\frac{DN}{1+\frac{NT_s}{z-1}},$$  then $$zY(z)+(NT_s-1)Y(z)=DN(zX(z)-X(z))$$
Use Z-transform formula $$z^{-n}Y(z) \leftrightarrow  y[k-n],$$ we have $$y[k+1]=DN(x[k+1]-x[k])-(NT_s-1)y[k]$$
Here x is the input and y is the output.
So my code for the PID differential part is
double y0=0.0;
double Ts=0.005;
double PID_DifferentialPart(double x,double D,double N){
  static double x_last,y_last;
  static int first=1;
  double y;
  if(first){
    first=0;
    x_last=x;
    y=y0;
  } else {
    y=D*N*(x-x_last)-(N*Ts-1.0)*y_last;
  }
  y_last=y;
  return y;
}

Using the same input (e.g. input sequence {0 2 3 4 2.5}, with D=2.0,N=100.0,Ts=0.005), the outputs are different. The output of my code is {0 400 800 1200 1100}. Output of Simulink code is {0 400 400 400 -100}. Where did I go wrong?


texts: Object Recognition using bag of features using matlab

texts: I want to write a MATLAB program for simple object recognition using bag of features. In short, I want to first extract the features from an image, create a visual library using those features, then cluster up the features belonging to one part together, hence creating different parts. now use these parts for matching. i have the basic idea but i don't know much math behind this, also i don't know how i can go about implementing it. can anyone help me with the matlab code for this?


texts: Understanding the $\mathcal Z$-transform

texts:
I was studying $\mathcal Z$-transforms and found pretty good material on the topic, though I feel I do not have a proper understanding of the concept. Could someone help me clarify this?

I know that we apply $\mathcal Z$-transforms to sequences and convert them to the $z$ domain and find the values of $z$ for which we get zeros and poles and also find the region of convergence. What I didn't get is:

Why are we multiplying some $z$ values to the given sequence and finding for what values of $z$ the result gets pretty huge or dampens out?
What does poles and zeros have to do with the given sequence? It's not like we are applying the transform to a system and finding the system's response, it's just a sequence.
And they also say that its a frequency domain transformation. How do we realise frequency here?


texts: How do I find the frequency of a particular bin in a cosine transform, and two other questions


texts:
I have this signal (as the text):
$$
x(t)=\sum_{k=-\infty}^\infty g(t-kT)\quad\text{with}\quad g(t)=\textrm{rect}(t)
$$
then we have an ideal bandpass filter $H(f)$ with frequency centered at $f=20\textrm{ kHz}$ and band $b=10\textrm{ kHz}$. Calculate the signal's output in frequency and in time. This is a periodic signal.
$$
X(f)=\sum_{k=-\infty}^\infty C_n\delta(f-nf_0)$$
doubts:

$$C_n=f_0\mathrm{sinc}(nf_0)\quad\text{and output}\quad\sum_{k=-25}^{-15}{C_n\delta(f-nf_0)}+\sum_{k=15}^{25}{C_n\delta(f-nf_0)}\quad ?$$


texts:
From Wikipedia:

  If the output signal slowly changes independent of the measured property, this is defined as drift

So my question is, if I have a gyroscope with sensor drift does that mean that $E\left\{f(x)\right\} \neq 0$ ? I have an excercise where it is written, that the sensor has drift but also that the $E\left\{f(x)\right\} = 0$. But to me one statement seems excludatory of the other one.


texts: Sparse recovery, Restricted Isometry Property for ILL-POSED problems

texts:
if $\mathbf x$ is $N\times 1$ sparse vector, and $\mathbf A$ is an $M\times N$ matrix with $M<<N$, and we measure $\mathbf y=\mathbf{Ax}$, then compressed sensing theory tells us that we can exactly recover $\mathbf x$ from $\mathbf y$ under some conditions for  $\mathbf A$ like RIP.

But what if $\mathbf A$ is an $M\times N$ matrix with $M>>N$ or $M=N$. Now the equation is not under-determined. But the equation can not be solved directly because $\mathbf A$ is now ill-conditioned matrix. To approximate to solution we can apply regularization. Since $\mathbf x$ still sparse we apply $l_1$ reguarization scheme as $$ ||y-Ax||_2^2 + \lambda ||x||_1$$ and find a sparse solution for x where sparsity level is dependent on $\lambda$.

The questions are

1) Do we still need to consider RIP for $\mathbf A$? What effects the acuuracy of the solution of $l_1$ regularization scheme except $\lambda$?

2) Why do we need RIP in compressive sensing and does it have any relation with the values of $M$ and$N$?

3) If $x$ is not sparse does this regularization fail for sure? Or it can still be work like $l_2$ regularization scheme. Maybe better?


texts: Sum of sines sound wave using microcontroller

texts:
I need to produce a signal which contains only several frequencies, for example, 17KHz, 17.2 kHz,17.6 kHz. The components I'm using are microcontroller (PIC16F1783 by microchip) which works at 80 kHz frequency, and produces an analog signal, which is then fed into an audio amplifier and then into a speaker.

First, I'm trying to output a single frequency, say 17 kHz. under timing and low memory constraints my algorithm is as follows:

I have an array of 400 elements which contains samples of a sine wave at constant intervals
Every clock cycle(80,000^-1 sec) I output sine_array[index] to the audio amplifier (and then to the speaker) where index starts at 0 and increased each cycle by: (17,000/80,000)*400=85.

However, for some reason, besides seeing the pure 17 kHz frequency, I also hear (and see on the oscilloscope) parasitic 5 kHz frequency!
I've also noticed that for frequencies such as 16 kHz (index increased by 80), 8 kHz there are no parasitic frequencies.
So I though it may have something to do with the fact that 400/80 is an integer, whereas 400/85 isn't.

Can anyone please suggest a creative way to sort this issue? The final goal is, as stated, to produce a signal which is composed of several frequencies with a rather small difference (around 200 Hz).


texts: Long Term Average Spectrum of large batch of audio files

texts: I'm trying to obtain the overall LTAS of a large batch of audio (wav) files. To do this, I should obtain the LTAS of each file individually, then take the average value across all LTASs for each frequency. I also need the standard deviation for each frequency. I'm doing this in Matlab, but am rather a noob. Can anyone offer any pointers to help me write a Matlab script to do this? Thanks!


texts: what is the relation between a minimum phase system and a FIR system?

texts: I want to see whats the relation between being an FIR system and a minimum phase system and when an FIR system can be minimum phase too with respect to have linear phase


texts: Highest rate achievable by current DSP and A/D chips?

texts:
I'm trying to determine whether an application I have in mind is feasible with currently available DSP technology (DSP chips and ADCs).

I want to eliminate an interfering signal using an LMS adaptive filter.  I have access to the source of the interference.

The bandwidth of the interfering signal is in the order of 1MHz  (maybe a little less, but just as a conservative assumption, let's say 1MHz).  The bandwidth of the signal I want to measure can be higher; possibly in the order of 10MHz or 20MHz.

EDIT: I estimate that the application will require sampling at no less than 12-bit resolution.

Does the above sound feasible?  Could you suggest some DSP chips that you'd estimate could potentially handle the above workload?

If it sounds like it's not feasible, can you give an estimate of what you think would be the highest bandwidth that current DSP chips could handle?

You may assume that the DSP chip will be doing exclusively this;  or maybe just some "lightweight" additional work.

Thanks!


texts:
Consider following LTI system $$y[n] - 2y[n-1] = x[n]$$ where $x[n]$ is the input to the system and $y[n]$ is the output. Let $x[n] = \cos[n\frac{\pi}{3}] + 2\cos[n\frac{\pi}{2} + \frac{\pi}{4}]$, determine the Fourier series coefficients for $y[n]$.

My try: First of all we should determine impulse response $h[n]$. So we have $h[n] - 2h[n-1] = \delta[n]$ and assuming initial rest condition leads to $h[n] = 2^nu[n]$. Then the frequency response could be computed as follows:
$$H(e^{j\omega}) = \sum_{n = -\infty}^{+\infty}h[n]e^{-jwn} = \sum_{n = 0}^{+\infty}2^ne^{-jwn} = \sum_{n = 0}^{+\infty}(2e^{-jw})^n$$
It's a geometric series and divergent because $|2e^{-j\omega}| = 2 >1$. So we should conclude that response to $x[n]$ doesn't exist? Also is initial rest a valid assumption here?

Edit: Let $a_k$ be the Fourier series coefficients for the input signal $x[n]$. Then the Fourier series coefficients for output is $b_k = a_kH(jk\omega_0)$ where $\omega_0$ is the Fundamental frequency of the input. The problem is that for each $k$ the geometric series diverges so $b_k$ doesn't exist.


texts: Why is interpolation a time varying system

texts:
I was reading about interpolation (Interpolation and Decimation of Digital Signals - A tutorial Review, Ronald E. Crochiere) and found that Interpolation filter is a time varying system. Can someone give an simple example why is it a time varying system. Because, as far as I understand, interpolation just stuffs zeroes between the samples to increase the sample rate by appropriate amount, even if we apply the input at a later time it should give the same output.
Also, even though the process is time varying we can represent it by using transfer functions. For instance, Zero order hold is represented by sinc filter. Why is this possible?


texts:
I was watching a YouTube tutorial on wavelets https://www.youtube.com/watch?v=wgRgodvU_Ms, where the presenter showed that if we multiply a Gaussian with a sine wave (and follow some other conditions), we get an example of a Morlet wavelet.

In the case of Haar, the scaling function (left) and the wavelet (right) have a different shape. How does one design a corresponding scaling function once a wavelet is designed? For example, a Mexican hat wavelet is a second derivative of a Gaussian with its signs reversed. How does the corresponding scaling function look like.?

A pointer to an introductory reference would be appreciated as well.


texts: Is a change in colorspace needed when comparing colors only evaluated by the computer?

texts:
I am confused about the need to change color space for color comparison. I have read about delta E, the Lab format, and I do understand that comparisons in the RGB color space will not seem appropriate to the human eye. However, my program uses a linear color scale to calculate velocity, from a color flow Doppler signal. It takes the mean color of a sample region and compares it to the colors of the scale to find its nearest neighbor using Euclidian distance. I do that entirely in the BGR (OpenCV) color space, as the example image below:

Here, I obtain seemingly correct velocity values for each color circle, but is it only by chance, or is my assumption correct that since the color comparisons take place internally, it does not matter what color space I am in?


texts: Time domain interpolation using FFT with zero padding on the end

texts:
I've got a situation where I'd like to use an FFT to do interpolation in time on some complex data (I need to go to the frequency domain anyways to window my data).

The notional way of doing this is to take an FFT, which, in the case of something like FFTW will put the frequency components "in-order" which can be viewed as having positive frequencies first, followed by negative.

$\mathcal{F}([\begin{array}{ccccccc} a & b & c & d & e & f & g & h \end{array}]) = \begin{array}{cccccccc} [A & B & C & D & E & F & G & H] \end{array}$

Where ABCD are positive frequency components, and EFGH are negative frequency components.

You can split the frequency response in two and zero-pad between D and E, but you have to be careful if your transform is an odd length, as well as possibly splitting the nyquist bin, and it's more cumbersome from a buffer-juggling perspective too.  It'd be better if I could pre-multiply by Fs/2, which will rotate the spectrum:

$\mathcal{F}(exp(j2\pi fs/2t) * [\begin{array}{ccccccc} a & b & c & d & e & f & g & h \end{array}]) = \begin{array}{cccccccc} [E & F & G & H & A & B & C & D] \end{array}$

Then it's just a matter of zero padding on the end, and getting back to the time domain.  This might be as simple as just transforming into a sufficiently large zeroed buffer.

But, it's the getting back to the time domain that's got me wrapped around the axle.  You can certainly inverse-FFT, and then multiply by -Fs/4 to shift things back to the way they should be, but I'd prefer to avoid the post-multiplication if I can (performance reasons).

My gut says that the shifted spectrum is now in the proper order (again in FFTW land) to just do a forward-FFT on it, but that might time-reverse the data, so perhaps a conjugate before transforming would be in order, but I'm not sure.  Can anyone help me puzzle this out?


texts: Deriving the Resolution Equation of an FMCW RADAR

texts:
I am interested in the equation that depicts the resolution of an FMCW radar. I write the same below.

$\Delta d = \frac{c}{2b}$

First I would like to have an intuitive explanation to the above derivation.

Secondly, if possible, I would like to know why the chirp's gradient doesn't play a role in determining the resolution of the FMCW radar. Does it not matter how fast I sweep the $BW$? If not so, what impact does the chirp gradient have?


texts: How condition for existence of Fourier transform is valid?

texts:
The condition for Discrete time Fourier transform to exist for function $f(n)$ is given as

$$\sum_{-\infty}^\infty |f(n)|  < \infty.$$

In case of continuous Fourier transform the difference is summation is replaced by integration in the above equation

It looks something
strange for me because on left side we are
taking the signal within 1. infinity
limits. Also,2. it is again of summation type and we are expecting right side to be less than infinity. How is it possible?

Let us take two  examples

  Suppose there is unit step signal $u(n)$
  with constant unity magnitude
  ranging from o to $\infty$ .then we
  can clearly see by putting value
  in the above formula

L.H.S=1+1+1+......upto ${\infty}$
So R.H.S. =$\infty$

So can I say that DTFT of unit step signal doesn't exist?

  Similarly, if I represent the  power signal like $cos(wn)$  using the exponential terms and compute it's value from the formula of condition of existence ,it will come out to be
  $\infty$ .

So can I say that DTFT of power signal like
 $cos(wn)$ doesn't exist?


texts: How can I perform timing recovery from the envelope of PSK signals?

texts:
I'm trying to reimplement FDMDV in GNU Radio, largely for educational purposes. It consists of 15 PSK subcarriers, 75 Hz apart, 50 symbols/sec each. The center one is BPSK and alternates phase each symbol, and the rest are QPSK and carry the data.

I'm somewhat stuck on how to perform timing recovery if this signal. Since each subcarrier has the same clock, it seems this should be worth something. The existing implementation seems to work by (full implementation at the end):

sum the envelopes of each subcarrier
this sum has a frequency component at the symbol rate, so determine the phase of that
drive this phase to a fixed value through a control loop

How might I implement this in GNU Radio? Is this something that can be accomplished with the standard blocks or will I need to implement a new block of my own? Or is there another method I could use to perform the same task?

function [rx_symbols rx_timing_M env fdmdv] = rx_est_timing(fdmdv, rx_filt, nin)
  samples_per_symbol = fdmdv.M;
  Nt = fdmdv.Nt;
  num_carriers = fdmdv.Nc;
  rx_filter_mem_timing = fdmdv.rx_filter_mem_timing;
  P = fdmdv.P;
  Nfilter = fdmdv.Nfilter;
  Nfiltertiming = fdmdv.Nfiltertiming;

  % nin  adjust
  % --------------------------------
  % 120  -1 (one less rate P sample)
  % 160   0 (nominal)
  % 200   1 (one more rate P sample)

  adjust = P - nin*P/samples_per_symbol;

  % update buffer of Nt rate P filtered symbols

  rx_filter_mem_timing(:,1:(Nt-1)*P+adjust) = rx_filter_mem_timing(:,P+1-adjust:Nt*P);
  rx_filter_mem_timing(:,(Nt-1)*P+1+adjust:Nt*P) = rx_filt(:,:);

  % sum envelopes of all carriers

  env = sum(abs(rx_filter_mem_timing(:,:))); % use all num_carriers+1 carriers for timing
  %env = abs(rx_filter_mem_timing(num_carriers+1,:));  % just use BPSK pilot
  [n m] = size(env);

  % The envelope has a frequency component at the symbol rate.  The
  % phase of this frequency component indicates the timing.  So work out
  % single DFT at frequency 2*pi/P

  x = env * exp(-j*2*pi*(0:m-1)/P)';

  norm_rx_timing = angle(x)/(2*pi);
  rx_timing = norm_rx_timing*P + P/4;
  if (rx_timing > P)
     rx_timing -= P;
  end
  if (rx_timing < -P)
     rx_timing += P;
  end

  % rx_filter_mem_timing contains Nt*P samples (Nt symbols at rate P),
  % where Nt is odd.  Lets use linear interpolation to resample in the
  % centre of the timing estimation window

  rx_timing += floor(Nt/2)*P;
  low_sample = floor(rx_timing);
  fract = rx_timing - low_sample;
  high_sample = ceil(rx_timing);
  %printf("rx_timing: %f low_sample: %f high_sample: %f fract: %f\n", rx_timing, low_sample, high_sample, fract);

  rx_symbols = rx_filter_mem_timing(:,low_sample)*(1-fract) + rx_filter_mem_timing(:,high_sample)*fract;
  % rx_symbols = rx_filter_mem_timing(:,high_sample+1);

  rx_timing_M = norm_rx_timing*samples_per_symbol;

  fdmdv.rx_filter_mem_timing = rx_filter_mem_timing;
endfunction


texts:
The signal of interest to me is the ocean wave height as a function of time at a particular location. Its power or variance spectrum is estimated by buoys at regular time intervals.

Specifically, the spectral data is discrete in both time and frequency. The spectral data spans 0.02 and 0.485 Hz, in approximately 0.02 Hz intervals. This spectrum is reported every hour (3600 seconds).

I would like to estimate the spectrum at finer time intervals (every 10 minutes, say). If have the spectrum at 1 am, 2 am, ... 11 am, what is the best way to estimate the spectrum at 5:30 am? I'm considering a linear interpolator, or a conventional decimator, but I'm wondering if there is an optimal way to do this.


texts:
I have implemented a texture recognition algorithm in MATLAB and it works well. Now, my next step is to summarize with the accuracy percentage. Now, I would like to summarize my work and attribute an accuracy to my algorithm.

I have around 1000 images and only one training image. The feature vectors are available for all the images in the dataset. Right now I am getting the top 5 matches from the dataset for a given image. Everything works well!

now i want to calculate accuracy for all the images in the dataset and average it out to get the final accuracy. Can someone tell me what could be the best method for calculating the accuracy and precision and how should i proceed in case of large datasets?


texts:
I am looking at data like this:

Source data: http://pastebin.com/raw.php?i=L6cd8d5K
The data is quasi-periodic, there is a distinct cycle but it is not an exact repeat, and the period changes considerably.  I would like to find accurate peak positions, and period durations, for each cycle.  Ideally the peak positions would match the intuitive notion of where the exact peak is.
What makes this hard:

noise: although it may not be entirely noise!  the variation is much greater when the values are high, and so the highest spikes tend to occur near where the "real" peak is.

the waveform is not symmetrical: neither left-to-right nor up-down.  doing something like a highpass filter will move peaks up and down by an amount that depends on their duration and the shape in between.

What is unlikely to work well:

bandpass filter combined with conventional peak picking: the position of the peaks will shift since values from the (asymmetrical) up and down slopes is added to where the peak is.  The really high values will not have a high weight in the result since the lowpass part of the bandpass filter will spread their energy to adjacent values.

autocorrelation: the shape of successive cycles is too different.

period durations measured from the leading (or trailing) edge of one cycle to the next: again shapes are too different.


texts: dealing with Huge Length (Order) FIR filters

texts:
I'm studying a project.

My signal is audio 8KHz Mono, 8 bits in Low frequencies 70Hz until 1200Hz.

When I perform the calculus with the transition band is very narrow (Transition Frequency) my Order results are (according to proposed window):

For this Example the worst case:

DeltaFrequency = 1.060180596543088 Hz
dNormalizedDeltaFreq= 1.060180596543088/8000

How is calculated the Order? My Answer is below, check the PD!

Rectangular: 6793
Hamming:24903
Hanning:23393
Kaiser:27353
Blackman:41503

If I try to use 16KHz, before orders will be doubled!

Then, I want to use Hardware only to perform the filtering FIR process, and I'm worry about of Time needed to perform this functionalities.

COMMENT1:

Reduce Sampling Frequency is not optional, I need to deal with 1200Hz signals too!

PD:  (This is not question!)

    if (windowType.toLowerCase().equals("rectangular")) {
      Size = (int)Math.ceil(0.9/dNormalizedDeltaFreq);
    }
    if (windowType.toLowerCase().equals("hamming")) {
      Size = (int)Math.ceil(3.3/dNormalizedDeltaFreq);
    }
    if (windowType.toLowerCase().equals("hanning")) {
      Size = (int)Math.ceil(3.1/dNormalizedDeltaFreq);
    }
    if (windowType.toLowerCase().equals("kaiser")) {
//      Mitra's book Handbook for Digital Signal Processing quotes Kaiser with a
//      simple estimate: N = [-20*log10(sqrt(dp*ds)) - 13]/[14.6(ws - wp)/2pi]
       double Att = 60.0;
      Size = (int)Math.ceil((Att - 7.95)/(14.36*dNormalizedDeltaFreq));
    }
    if (windowType.toLowerCase().equals("blackman")) {
      Size = (int)Math.ceil(5.5/dNormalizedDeltaFreq);
    }

QUESTION 1: it is feasible to use a FIR filter of this Order?

QUESTION 2: Which single device (hardware) help me perform this function?

QUESTION 3: How could be calculated the delay time?


texts:
I am new in matlab and signal processing. The time series that have been used are obtained from accelerometer in a building. As far as I understand both the time series' length and window function length(hamming or hanning) should be same. But for smoothing the fourier amplitude spectrum, optimum window length can vary. How this procedure should be apply in matlab?

How can apply a window function that is shorter than my time series.

Thanks for reading!


texts: Estimate Image power spectrum of the noise

texts:
I want to implement Wiener filter. I have original image and image blured with linear motion. The one thing I am lack of(to implement it) is Noise power spectrum.

I know that it is hard task. I've took a look on different papers http://kmh-lanl.hansonhub.com/publications/medim98a.pdf and https://www.ima.umn.edu/2005-2006/MM8.9-18.06/activities/Wohlberg-Brendt/team5_rep.pdf

I don't understand why it is that complicated to just find noise spectrum if I have original signal(image)?

If I take noise image minus original image, it will be my error image. And if I estimate variance of it, that will be my Noise variance(Noise Power)?

I totally confused and here is my question:
" How to estimate Power spectrum of the noise of the blured image if I have original image?"

Btw, I don't want get some matlab function which will do all of it, but rather algorithm or equation that I can implement.

Happy New Year!


texts: FFT of square wave - what does output represent?

texts:
I am really new to FFT and signal processing. I am doing an analysis of square waves with FFT and I am trying to understand why the FFT output on the frequency domain has a downward slope for square waves. Specifically:

Why in FFT output for a square wave, bin 0 has the highest amplitude and bin 1 has the second highest amplitudes no matter what the original frequency of the input signal in the time series domain representation is.
I read that bin 0 has to do with the ADC and that is the sum of the whole signal and that is why it is 0Hz. What does that mean?
What does amplitude represent in FFT output like below?


texts: Exponential decaying step response in LTI System

texts:
I'm attempting to better understand the relationship between step responses, impulse responses, and convolutions. Say that I have a system where if I apply a constant input, my output decays from a certain value to an equilibrium value. This gives me the step response of a system that I believe to be LTI based on the system's differential equation.

PDE:

$$\frac{\partial^2u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}-\frac{u}{r^2}=\frac{1}{H_A k}\frac{\partial u}{\partial t}+\frac{1}{H_A k}\frac{r}{2}\frac{\partial \epsilon}{\partial t} $$

ODE in S Domain

$$\frac{\partial ^2 \bar{u}}{\partial \hat{r}^2}+ \frac{1}{\hat{r}}\frac{\partial \bar{u}}{\partial \hat{r}}-\frac{\bar{u}}{\hat{r}^2}-s\bar{u}=\frac{s\bar{\epsilon}\hat{r}}{2} $$

My specific step response is fairly complex and I went more in-depth into this specific system in another post, but a simpler version that still captures my issue could be given as:

$$
Y(t) = 1+e^{-0.1t}
$$

In an LTI System, it's my understanding that the derivative of the step response yields the impulse response, which can then be used to obtain the output of the system for any arbitrary input through convolution.

$$
h(t) = \frac{dY(t)}{dt} = -0.1e^{-0.1t}
$$

So now, if I have this impulse response, which is negative, and I attempt to convolve this with a unit step function to recover my original Y(t), I get something that is both the wrong shape and the wrong sign. Can anyone help point out my mistake?


texts: Closed form of $\mathcal Z$-transform : decomposition signal $x(n)$

texts:
The text of my exercise ask :

  Determine the closed form of the $\mathcal{Z}$-transform for this $x(n)$  $$ x(n) = \begin{cases} |n-N|  & \text{if 0<$n$<2N} \\ 0 & \text{elsewhere} \end{cases} $$

I see $x(n)$ like a composition of two signals

$$ x(n)= x'(n)+x''(n) $$
I use for greater semplicity $x_1(n)$
$$x_1(n)= (u(n)-(u-(N+1))$$
whence
$$x'(n)= (n-N) x_1(n) $$
$$x''(n)= x'(n-N)$$

Determination $X(Z)$ it easy, but if
the above equations are correct, for this reason i need post this question here.

$$X_1(Z)= \frac{z^{-(N+1)}-1}{Z^{-1} -1}$$
$$X'(Z)= z^{-1} \frac{{d}}{dz^{-1}} X_1(Z)- N X_1(Z)$$
$$X''(Z)= X_1 (\frac{1}{Z})z^{-N}$$

This is my solution.

$$X(Z)= X'(Z)+X''(Z)$$


texts: State Space model of a differential eqn. for use in Least Squares

texts:
The code in this page solves the Least Squares problem for the following dynamic model:

$\dot{y}=ay+bu$

where $a$ and $b$ are constants, $u$ is an input. The code is as follow:

t=[0:0.1:10];
m=length(t);
y=zeros(m,1);
a=-1;b=1;c=1;d=0;u=[100;zeros(m-1,1)]; % u: impulse with magnitude of 100
[ad,bd]=c2d(a,b,0.1);
y=dlsim(ad,bd,c,d,u);
ym=y+0.08*randn(m,1); %adding Gaussian noise

w=inv(0.08^2);
h=[ym(1:m-1) u(1:m-1)];
xe=inv(h'*w*h)*h'*w*ym(2:m)

ye=dlsim(xe(1),xe(2),c,d,u);

What's is not clear to me is why $a=-1$ and $b=1$. The model as given does not specify a value for either of $a$ and $b$. ($c=1$ and $d=0$ are fine, however).


texts:
I have a C# desktop application. I am using EMGU as a wrapper for OpenCV.

I am receiving images in JPEG format that were originally taken by an analogue camera.

The quality of the image is poor (compared to digital cameras).

I want to see if I can use filters/methods derived from OpenCV/EMGU that could enhance the image.

What I mean by enhancement is to make the image a bit sharper.

I have played around with median, bilateral, contrast stretching to see if the image can be improved.

Nothing has worked.

Even if it can be confirmed to me that nothing I can do will enhance the image that will be an answer in itself. I can stop trying and persuade the client to use a digital camera instead...

This is the image:

This is the code that I apply to the original image to get the image you can see here:

        int factor = Shared.factor;
        int xVal = 288;
        int yVal = 360;
        using (Image<Bgr, Byte> color = new Image<Bgr, byte>(grid))
        {
            byte[] data = null;

              for (int x = 0; x < xVal; x++)
              {
                  for (int y = 0; y < yVal; y++)
                  {
                      grid[x, y, 0] = (byte)(grid[x, y, 0] / factor);
                      grid[x, y, 1] = (byte)(grid[x, y, 1] / factor);
                      grid[x, y, 2] = (byte)(grid[x, y, 2] / factor);
                  }
                }
                color.Data = grid;
                data = ImageToByteArray(color.ToBitmap());
        }

    private byte[] ImageToByteArray(Image imageIn)
    {
        byte[] data = null;
        try
        {
            using (MemoryStream ms = new MemoryStream())
            {
                imageIn.Save(ms, System.Drawing.Imaging.ImageFormat.Jpeg);
                imageIn.Dispose();
                data = ms.ToArray();
            }
        }
        catch (Exception ex)
        {
         //do something
        }
        finally
        {
            if (imageIn != null)
            {
                imageIn.Dispose();
            }
        }
        return data;
    }

As requested this is an original image but from a different view:


texts:
Reading this paper about WSOLA (Verhelst, Roelands) (PDF)

At page 2 there's a "time-warped STFT's synthesis function" involving $v(n)=w^2(n)$ which is a windowing function. A windowing function squared.

What's the interpretation/meaning of windowing function squared? What does it mean visually?


texts:
We have designed a system which has an 8X8 array of photo diodes, we sense 10ns pulses of laser falling on those photo diodes with a repetetion rate of maximum 50Khz i.e. 20us, we were successfully able to sense them and detect them,please note that our system only senses those laser pulses and gives out the frequency of those pulses and amplitude/energy of those pulses, but doesn't do any kind of imaging

but now the problem statement has become a bit complex, for increased requirements we are forced to switch to a 64X64 array model, now we have to have 4096 stages of processing which will make the design very bulky

now i am thinking of solving this problem with a different approach, so in this need i need some suggestions, what alternatives can solve the problem by reducing the electronics stages

i thought of solving it by using a line scanning mechanism where 64 stages of processing will be given to 64 rows and 64 columns will be switched to use the same stages, but with this it will be very difficult to sense the 10ns pulses, because my sampling rate for each row should be at least 2GHz to acquire at least 10samples of the pulse the total sampling rate for scanning all columns is 2GHz X 64 = 128GHz sampling rate, working at this rate is very costly affair and over ambitious to go at a frame rate of 2GHz

please throw some light by suggesting any alternate techniques that these to solve the problem of increased replica of electronic stages


texts: How to apply baseline methods to MNE preprocessed raw EEG signals?

texts:
I am new at EEG preprocessing and classification. I have followed Clemens Brunner's blog https://cbrnr.github.io/2018/01/29/removing-eog-ica/ as a tutorial for MNE.

This is how I have preprocessed and singled out all the individual trials for all subjects in the dataset BCI Competition IV Dataset 2a. The dataset can be found here: http://bnci-horizon-2020.eu/database/data-sets

ch_names = ['Fz', 'FC3', 'FC1', 'FCz', 'FC2', 'FC4', 'C5', 'C3', 'C1', 'Cz',
        'C2', 'C4', 'C6', 'CP3', 'CP1', 'CPz', 'CP2', 'CP4', 'P1', 'Pz',
        'P2', 'POz']

mat = loadmat('A01T.mat')

for j in range(6):

    timelist = mat['data'][0][j+3][0][0][1]    //this gives starting of each event
    timelist = list(timelist)
    timelist.append(np.array([96735], dtype='int32'))    //length of eeg for single channel is added
    timelist = np.array(timelist)

    actionlist = mat['data'][0][j+3][0][0][2]    //these are the labels

    eeg = ((mat["data"][0][j+3][0][0][0]*10e-6).T)[:22]
    raw = mne.io.RawArray(eeg, info)
    raw.set_montage("standard_1020")
    raw.filter(1, 30)

    raw_temp = raw.copy()
    ica = mne.preprocessing.ICA(method='infomax', fit_params=dict(extended=True))
    ica = ica.fit(raw_temp, picks=['eeg'])
    ica.detect_artifacts(raw_temp)
    raw = ica.apply(raw)
    trial = raw.get_data()

    label = []
    train = []

    for k in range(48):
        trial_temp = trial[:, timelist[k][0]:timelist[k+1][0]]
        train.append(trial_temp)
        label.append(actionlist[k][0] - 1)

The file is saved in the shape of 22 channels and corresponding measures for the next 7-8 seconds. Frequency is 250Hz.

How do I apply baseline methods such as SVM to this preprocessed data?

Also, how do I calculate the power of the signals using FFT or any other method?


texts:
I wrote a python program to generate OFDM packets (i.e. the IQ samples) with a custom payload to, for example, be transmitted via an SDR.
After mapping the chunk of bits corresponding to one ofdm symbol to a constellation mapping (e.g. 64QAM), I perform the IFFT (using numpy.fft.ifft) and add the 16-sample cyclic prefix.
In the IEEE standard, the resulting 80-sample time snippet is multiplied by a windowing function (x0.5 for first and last value). However, when I test the whole thing and perform FFT on the 64-sample snipped again the whole mapping has a rather large error.

If I omit the windowing its perfect. I can't really understand the heavy influence, since its only one single value that has been altered by the window function in the end (concerning the fft).

Now I've read in the standard that symbols are appended one after another with one sample overlap, which I am currently not doing, so I guess there lies my mistake. If I append them with an overlap of one though, I would end up with a wrong overal sample count, wouldn't i? ( i.e. $80 * nSymbols - nSymbols$ instead of $80 * nSymbols$)

Thanks for every hint!
Cheers Lucas


texts: How to calculate the plane at infinity? either in Matlab or Python

texts:
I tried to implement this function to compute the plane at infinity $n_\Pi$ from the equation:

where $w^{-1}$ is the image of the absolute conic IAC, $M_j$ is the Projection Matrix of image $j$,
but I'm not sure if this implementation is right, because I don't have a good idea about symbolic expressions and how they are solved.
so I just mimicked other codes, but I'm still confused, because when I used this part to calculate the intrinsic matrix, the results of the intrinsic matrix weren't converging.
can you please tell me if there are issues in my implementation, and how to fix it please?

The whole code, and data is here

Thanks in Advance

function P_infty = PlaneAtInfinity(K,P)

% Symbolic variables

X =  sym('X', 'real');
Y =  sym('Y', 'real');
Z =  sym('Z', 'real');
L2 = sym('L2','real');

N = [X; Y; Z];
W_inv = K*K';
% Quadric
Q = [W_inv       ,   (W_inv * N)    ;
    (N' * W_inv) , (N' * W_inv * N)];

% Autocalibration equation
Calibration = P * Q * P';
% func = cross(W_inv, Calibration);
% solution = Optimize(zeros(3,3),@func);
% solve linear equations
solution = solve(Calibration(1, 1) == (L2 * W_inv(1, 1)), ...
    Calibration(2, 2) == (L2 * W_inv(2, 2)), ...
    Calibration(3, 3) == (L2 * W_inv(3, 3)), ...
    Calibration(1, 3) == (L2 * W_inv(1, 3)));

P_infty = [double(solution.X(1));...
    double(solution.Y(1));...
    double(solution.Z(1))];
end

and the main script:

P = [4.35310807409356, 1.44651724424134, 36.0849770394682, 6.6274893804074;
22.9809618064503, 8.56379922408173, -2.94837061198918, 36.6987313462383;
-0.00438517470946907, -0.00160149495653558, 0.241477646472821, -0.00696180145320084];

K = [870.258271610840, 0, 279.100746891319;
0, 812.178628533806, 260.930050222247;
0, 0, 1];
PlaneAtInfinity(K,P)


texts: Create Matched Filter Reference Signal when Scaling Unknown

texts:
I am working on a problem which raised what should be a simple question, but I'm stuck on it. The objective is to detect a shifted and scaled version of a known signal in the presence of Additive White Gaussian Noise (AWGN).

Think of a radar system. In a radar system, you know the transmitted signal. The returned system is delayed in time and scaled. I understand that a time-delay, only shifts the location of the matched filter's peak output. No problems there. But if the scaling is unknown, then how does one construct the reference signal?

Using a reference signal with the an amplitude that doesn't match the amplitude of the signal in the noisy data will mean that a threshold tests based on that reference signal will be incorrect.

I'm guessing that in a radar system, the noise power is approximately known, and so the power of the received signal is estimated by taking the power of the received signal and subtracting the noise power. That value could then be used to scale the reference signal. Still, that would only be and approximation, and it would still effect the probability of detection and false alarm in some way. Any help is appreciated. Thanks!


texts:
What are the different constraints for mastering audio content for either Vinyl or CD?

By the way, I'd like to know if there're any papers or documents, with a comparison between  different digital format (mp3, wav, ...) vs vinyl.

Thank you.


texts: Why the plots of spline and cubic interpolation are exactly same?

texts:
I am trying to watch difference between cubic interpolation and spline interpolation using matlab plot but i am getting same plots in both cases using interp1 command
My code is below. In below code ,if i use spline in place of cubic, i still exactly get same plot,why?
clc
clear all
close all
x = 0:pi/4:2*pi;
v = sin(x);
xq = 0:pi/16:2*pi;
vq2 = interp1(x,v,xq,'cubic');
plot(x,v,'o',xq,vq2,':.');


texts: Unit Problem in Designing a Filter for a Given Auto Correlation Function

texts:
Given a WSS process with the following Auto Correlation Function:

$$ r\left ( \tau \right ) = {\sigma}^{2} {e}^{-\alpha \left | \tau \right |} $$

The Laplace Transform would be:

$$ R \left ( s \right ) = \mathfrak{L} \left \{ r \left ( \tau \right ) \right \} = \frac{-2 \alpha {\sigma} ^ {2}}{\left ( s - \alpha \right ) \left ( s + \alpha \right )} $$

Hence the filter would be of the form:

$$ H \left( s \right) = \frac{c}{s + \alpha} $$

My question is about units.
In the Auto Correlation Function the units of $ \alpha $ are [Hz].
While in the filter form, assuming $ s = j \omega $ the units are [Rad / Sec].
How this conflict can be resolved?
What a I missing?

Thanks.


texts: Can you give example applications using single bin sliding DFT?

texts: I need to detect a single tone signal in real-time. The amplitude and the phase of the single tone signal are changing in time. I want to detect both the amplitude and the phase of the signal. I do not want to do FFT since I only need to see a single frequency bin and also FFT is not suited well for real-time applications. I see that there are sliding Goertzel and sliding single bin DFT methods. These methods seem easy to apply and both need very fewer computations compared to DFT. However, they seem like have stability issues. Some damping factors are used to solve the stability problems. However, I still have doubts on the sliding DFT method. Is there any point that I need to pay attention while applying single bin sliding DFT? Can you give examples of real world applications that are using single bin sliding DFT?


texts:
Let $T(s)$ be a transfer function that describes a mechanical system, where the input is force and the output is position.
And let $[A,B,C,D]$ be the equivalent state-space representation of $T(s)$, where:
$$
    \dot{x}(t) = Ax(t) + Bu(t) \\
    y(t) = Cx(t) + Du(t)
$$
and let $[A_d,B_d,C_d,D_d]$ be the discretized model of $T(s)$, using zero-order hold on the inputs:
$$
    x_d(k+1) = A_dx_d(k) + B_du(k) \\
    y(k) = C_dx_d(k) + D_du(k)
$$
I know there are infinite possible realizations of $[A,B,C,D]$ that represent the same transfer function $T(s)$.

How can I force the state vector, $x(t)$, to have a simple physical meaning? in this example, since it is a mechanical system, the elements of the state vector $x(t)$ would probably be position, velocity and acceleration. I can guess the matrix $A$ must have this form:
$$
A = \left[\begin{array}{ccc} * & * & * \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{array}\right]
$$
to express the derivative relationship between the state vector elements.
Same question, but for the discrete time case: How can I force the state vector, $x_d(k)$, to have a simple physical meaning? In this case, I'm not sure how "velocity" and "acceleration" can be expressed.


texts: Image classification using SIFT features and SVM

texts:
I am hoping someone can explain how to use the bag of words model to perform image classification using SIFT/SURF/ORB features and a support vector machine?

At the moment I can compute the SIFT feature vectors for an image, and have implemented a SVM, however am finding it hard to understand the literature on how use the bag of words model to 'vector quantize' the SIFT features and build histograms that give fixed size vectors, that can be used to train and test the SVM.

Any links to tutorials or literature on the topic is welcome, thanks


texts: Confusion in the derivation of convolution formula

texts:
$$x(t)=\int_{-\infty}^{\infty}x(\tau)\delta(\tau-t)d\tau=\int_{-\infty}^{\infty}x(\tau)\delta(t-\tau)d\tau$$
Due to linearity,
$$y(t)=T\{x(t)\}=\int_{-\infty}^{\infty}x(\tau)T\{\delta(\tau-t)\}d\tau=\int_{-\infty}^{\infty}x(\tau)T\{\delta(t-\tau)\}d\tau$$
considering the impulse response of system to be h(t),due to time invariance of unit impulse signal ,
$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(\tau-t)d\tau=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau$$
here, i know i have done something wrong, because
$y(t)=\int_{-\infty}^{\infty}x(\tau)h(\tau-t)d\tau$
is the formula for cross-correlation between x(t) and h(t)
and $y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau$ is the formula for convolution between x(t) and h(t)
but i cannot find, what exactly i have done wrong. could you help?


texts:
I am applying 8 samples and triangular window to my fft engine. I am getting the results but not accurate (approximately equal with some tolerance).

I am calculating MATLAB reference as below :

a=[1,2,3,4,5,6,7,8];
b=triang(8);
c=a*b;
d=fft(c)

Is it the correct way to apply window for signal processing unit?


texts: Intuition behind impulse response terms in convolution

texts:
We're learning about convolution in my signals and systems class right now. I have been able to do all of the problems by simply working out the respective sum/integral, but I'm still having trouble gaining the intuition behind it.

Consider the following example. Let $x[n]$ be a discrete-time signal and input it into some LTI system with impulse response $h[n]$. Then,

$$
y[n] = \sum_{k=-\infty}^{\infty} x[k] h[n-k].
$$

Let's plug in some values to make this more concrete. Suppose that we want to compute $y[6]$. Well,

$$
y[6] = \sum_{k=-\infty}^{\infty} x[k] h[6-k] = \cdots + x[4]h[2] + x[5]h[1] + x[6]h[0] + x[7]h[-1] + \cdots.
$$

I understand the shifting, but I feel as if the multiplications should be in a different order. Namely, why are we multiplying $x[7]$ by $h[-1]$. I feel as if we should be multiplying it by $h[1]$, since we've essentially shifted everything to the right by $6$ units, to $6$ is the new $0$, which would mean that $7$ is the new $1$ (under the shifting). I have it backwards, and understand why when I work out the math, but why?

More concretely, my question is the following:

With respect to the above example, what exactly is the meaning of $h[2]$? $h[-1]$? $h[k]$ in general?


texts:
Wiener–Khinchin theorem states that autocorrelation function and power spectral density are a Fourier-transform pair - see Wikipidia (and lots of other resources).
That means autocorrelation should be able to be obtained by inverse Fourier transform the spectrum. The following code (run within Octave, with "pkg load signal") shows the Fourier transform of the autocorrelation DOES look like the spectrum, but the inverse Fourier transform of the spectrum does not look like the autocrrelation. What did I do wrong?
### A signal's autocorrelation and its Engergy Spectral Density are Fourier transform pairs.

###   signal
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t)+randn(size(t));
#plot(x);

###    autocorrelation
Rxx = xcorr(x);
#figure();
plot(Rxx); title("Rxx");

###    autocorrelation FT
RxxDftAbs = abs(fftshift(fft(Rxx)));
freq = -Fs/2:Fs/length(Rxx):Fs/2-(Fs/length(Rxx));
figure();
plot(freq,RxxDftAbs); title("RxxDftAbs");

###    Energy Spectral Density
xdft = abs(fftshift(fft(x)));
x_esd = xdft.^2;   # ESD is the same as autocorrelation FT. Here for visualization purpose, using absolute values.
freq = -Fs/2:Fs/length(x_esd):Fs/2-(Fs/length(x_esd));
figure();
plot(freq,x_esd); title("x esd");

### ?????????????
### is it possible to get autocorrelation from the ESD by inverse Fourier transform?
### ?????????????

###    IFT of ESD
x_esd_idft_abs = abs(ifft(fftshift(x_esd)));
figure();
plot(x_esd_idft_abs); title("x esd ift");

Thanks in advance.


texts: How Are Images Considered Non Stationary Signal When They Are Invariant to Time?

texts: I have read Wavelets are better than Fourier in dealing with non-stationary signals such as images, but I don't understand how images are considered stationary??


texts: How does MATLAB's $\tt upfirdn$ select extra points when downsampling?

texts:
I had a signal 512 samples long and I used upfirdn to filter and downsample it by a factor of 2. I thought this would result in a 256 sample long variable, but instead it comes out as 261.

The MATLAB help says this:

  The length of the output, yout, is
  ceil(((length(xin)-1)*p+length(h))/q)

(where xin is the signal, h is the filter, p is the upsample ratio, and q is the downsample ratio) which, if you work it out (filter length is 10), does in fact come out to 261.

What I want to know is how it selects points when downsampling so that I can reproduce the upfirdn behavior in another language (long story, not important). Where does that formula come from?


texts: A wide sense stationary random process that is not second order stationary

texts:
I have been reading Peebles Probability, Random Variables, and Random Signal Principles and it claims that second-order stationarity is sufficient to guarantee:

$E[X(t)]$ is a constant
$R_{XX}(t1,t2) = E[X(t2)X(t2)]$ depends only on the difference $t1-t2$ of the two arguments and not on the individual arguments $t1$ and $t2$.

However it says that second-order stationarity is often more restrictive than necessary and hence the concept of a process that just complies with this conditions, named wide-sense stationary process, was created.

A second-order stationary process is defined as having a second order density function that is only dependent of time difference.

I would like to know:

What is meant by "more restrictive than necessary"? Is it just having to prove that second order density function is dependent only of time difference? What other aspects does a second-order stationary process has that a wide-sense stationary process doesn't?
If a second-order stationary process is also first-order stationary. What would a wide-sense stationary process require to be first-order stationary? A first-order density function that does not change with a shift in time?
How high must the N in Nth-order stationary process be to be considered a Strict-Sense Stationary process? Infinite?


texts: Gradient echo MRI. Do the gradients have to be the same duration?

texts:
I currently got hold of the book "Handbook of MRI Pulse Sequences" and am looking into computing the duration, gradient strength and ramp times for the phase encoding gradient, readout prehasing and the slice prephasing gradients.

Now, every GE pulse sequence diagram that I look at, shows these gradients to be of the same duration i.e. they start and end at the same time during the sequence. Now, going through some of the calculations in the book for a scanner with some maximum gradient strength, slew rate, receiver BW, field of view and the k-space matrix size, I can compute the optimal gradient strengths and the durations for these gradients.

However, the optimal gradients are not the same duration. Does this mean I have made a mistake in my calculations or do they overlap these gradients in the pulse sequence diagrams for convenience?


texts: What Is the Best First Order IIR (AR Filter) Approximation to a Moving Average Filter (FIR Filter)?

texts:
Assume the following first order IIR Filter:

$$ y[n] = \alpha x[n] + (1 - \alpha) y[n - 1] $$

How can I choose the parameter $ \alpha $ s.t. the IIR approximates as good as possible the FIR which is the arithmetic mean of the last $ k $ samples:

$$ z[n] = \frac{1}{k}x[n] + \frac{1}{k}x[n-1] + \ldots + \frac{1}{k}x[n-k+1] $$

Where $ n \in [k, \infty) $, meaning the input for the IIR might be longer than $ k $ and yet I'd like to have the best approximation of the mean of last $ k $ inputs.

I know the IIR has infinite impulse response, hence I'm looking for the best approximation.
I'd be happy for analytic solution whether it is for $ {L}_{2} $ or $ {L}_{1} $ cost function.

How could this optimization problems can be solved given only 1st order IIR.

Thanks.


texts:
I have a gyroscope which reports x,y,z values to my embedded device.
The rate that one sample appears is variable and has a jitter of around 10ms.

However I later want to have a log file with a fixed sample rate (higher than the input sample rate). I found this post What is an algorithm to re-sample from a variable rate to a fixed rate?

Unfortunately it seems way too complicated compared to what I want to do. For me it would be fine simply interpolating linearly between the samples, but I'm missing the correct algorithm to do so. Is there any common way on how to interpolate a variable sample rate to a fixed sample rate?


texts: Extracting amplitudes of components from FFT

texts: I have a 4096-point FFT of a continuous pipe organ sound (no windowing).  I wish to extract the relative amplitudes of the various harmonically-related sinusoidal components for voicing analysis.  I can easily find the bins with the local maximum power for each component, but since the frequencies are not exactly in the middle of a bin, some of these peaks are spread between two adjacent bins, and the amplitude in each bin is less than it would be if the frequency were shifted to exactly hit that bin.  For frequency analysis I have used various interpolation formulas with great success.  But now I am doing amplitude analysis, and I want to know if there is a similar interpolation formula to estimate the amplitude of a component in a way that is independent of whether that component is exactly in the center of a frequency bin or is between two bins.  My current best guess is the take the sum of the powers (magnitude squared) of the bin with the peak power together with the bins immediately adjacent to it.  The sum of those three bin powers would then represent the power of that one component of the sound.  Is there a better way?


texts:
I am trying to filter a sinusoid ($f_{sample}=50 Hz$) ,, i need to remove all frequencies above 10 Hz with accepted transition region 2 Hz (stop band freq=12 Hz)

the question is why after filtering this signal is the 20 Hz still appearing in the FFT? is there something wrong ? I know results should be the same in frequency and time domain filtering but the first in multiplication and the 2nd in convolution.

here is the code :

 t=[0:1/50:1];
 y=cos(2*pi*20*t)+0.5*sin(2*pi*5*t);
 fsamp=50;
 frange=[10 12];
 mag=[1 0];
 devs=[0.02 0.02];
 frange=[10 12];
 [n,fo,ao,w] = firpmord(frange,mag,devs,fsamp);
 n = n + rem(n,2);
 b=firpm(n,fo,ao,w,'hilbert');
 %%% filter in time domain
 final1=filter(b,1,y);
 finalfft=abs(fft(final1));
 t=[0:1/50:(length(finalfft)-1)/50];
 f=[(t/t(end))*50];
 figure(1);plot(f,abs(fft(y)));
 figure(2);plot(f,finalfft);

 %%% filter in freq domain
 paddedfft=fft(y);
 bfft=fft(b,length(y));
 finalfft=paddedfft.*bfft;
 final1=ifft(finalfft);
 t=[0:1/50:(length(finalfft)-1)/50];
 f=[(t/t(end))*50];
 t=[0:1/50:(length(final1)-1)/50];
 figure(3);plot(f,abs(fft(final1)));


texts: Equal contribution criteria during image pyramid construction

texts:
In their paper on Image Pyramids, Burt and Adelson state the following about pyramid generation.

  An additional constraint is called equal contribution. This stipulates that all nodes at a particular level must contribute the same total weight($=\frac{1}{4}$) to nodes at the next higher level. Let $\hat{w}(0)=a$,$\hat{w}(-1)=\hat{w}(1)=b$,$\hat{w}(-2)=\hat{w}(2)=c$. In this case, equal contribution requires that $a+2c = 2b$.

How is this weight condition derived and why is the weight $\frac{1}{4}$? For a larger context of the question, the screenshot of the relevant section in the paper can be seen below


texts:
I am working on kind of signals that can be classified as 1/f Noise.
In this article is explained that the time history of 1/f noise follows a power law proportional to $t^{-\alpha}$. I am very surprised about this, because when I think to this kind of signal histories, they do not seems to me proportional to $t^{-\alpha}$ or a series of $t^{-\alpha}$.
Have someone an explanation for this?
I have also read this article, that highlights the proportionality to $t^{-\alpha}$ and introduces the theory of the "Self-Organized criticality". But I didn't find an explanation to my problem.


texts:
Given a continuous LTI system with transfer function
$$H(s)= -\frac{2s}{(s+6)(s+2)}$$

Plot the location of the pole(s) and zero(s)
Find all possible regions of convergence
From the problem above find the impulse response

Here's my attempt

The poles are $s_1=-6$ and $s_2=-2$ and the zero is $s_3=0$

The all possible roc-s are

Using inverse laplace transform I found $h(t)=e^{-2t}-3e^{-6t}$. For the right-sided signal I multiply $h(t)$ with $u(t)$ and I got $h_{RS}(t)=(e^{-2t}-3e^{-6t})u(t)$ and for the left-sided signal I multiply $h(t)$ with $u(-t)$ and I got $h_{LS}(t)=(e^{-2t}-3e^{-6t})u(-t)$.

My confusion is at the two-sided signal. My solution is $h_{TS}(t)=e^{-2t}-3e^{-6t}$ and my friend's solution is $h_{TS}(t)=-e^{-2t}u(-t) - 3e^{-6t}u(t)$. When I plotted the signal both are a two-sided signal. Since my professor didn't submit the solution we didn't know who's right or wrong.
Thank you.


texts: Correlating a sequence or a code with itself at QPSK symbol rate

texts: I am sending a QPSK signal that represents a certain sequence or code and  that needs to be correlated with itself on the receiver (there are other sequences/codes in the system, and they are all orthogonal). The result of the correlation should  indicate which code was sent . Suppose an AWGN channel, and at the output of the I and Q matched filters on the receiver  I have one sample/symbol that is a 'sufficient statistic" , so one sample on I and one sample on Q , and suppose timing synchronization is perfect and no other impairments besides AWGN (SNR is high enough), my question is the following:  Can correlation take place at the symbol rate ? or does one need more than one complex sample/ symbol in order to do the correlation?


texts: Computing cross-correlation of two images using OpenCV

texts:
I have a program that uses OpenCV to compute either the convolution or cross-correlation of an image with a specified kernel. I compute cross-correlation by setting the conjB flag to true when calling cv::mulSpectrums.

Convolution consistently produces correct results. However, the result of cross-correlation is always wrong.

Based on the convolution results, the spectra of the kernel and image are correct. The only change is that the kernel spectrum is conjugated before multiplying. Both were computed using cv::dft, and so are in the CCS-packed format expected by mulSpectrums.

I'm comparing the result with the output from scipy.signal.correlate, using Python. I tried flipping the the kernel before computing the DFT, and that produces the correct cross-correlation.

cv::dft( dft_buf, dft_buf, 1, img_in.rows );
cv::mulSpectrums( dft_buf, kernel_dft, dft_buf, 0, CORRELATE );
dft( dft_buf, dft_buf, cv::DFT_INVERSE + cv::DFT_SCALE, kernel.rows/2 + img_in.rows );

int col_offset = ( kernel.cols % 2 ) ? ( kernel.cols/2 ) : ( kernel.cols/2 - 1 ); // is odd?
int row_offset = ( kernel.rows % 2 ) ? ( kernel.rows/2 ) : ( kernel.rows/2 - 1 );

cv::Mat result_roi( dft_buf, cv::Rect( col_offset, row_offset, img_in.cols, img_in.rows ) );
result_roi.copyTo( img_out );


texts: How do I get the phase angle from Cross Correlation?

texts:
Im trying to cross correlate two signals in matlab and get the phase difference between the signals.

For cross correlation (the idea is to do it without xcorr) I used:

Cxx=fftshift(ifft(fft(x,N).*conj(fft(y,N))))/(norm(x) * norm(y));

I get the result and not sure about the reference point I have to take for phase calculation Now which point should I take as zero on time scale and how do I get the phase difference from this result.

Thanks.


texts: Extracting ridges in automatically in image

texts:
I am trying to extract ridges from an image, calculated using the eigenvalues of the Hessian matrix (see answer in Best way of segmenting veins in leaves?). In this image I have zeroed out the negative eigenvalues.

Visually I can see the regions I want, but how would I be able to sensibly extract them automatically? If I LogNorm() the colorbar it is not so clear

Initially I thought if the noise in the image can be estimated, I can mask the image accordingly, but the only estimate to the noise I can find is its standard deviation (assuming Gaussian noise) in skimage.restoration.estimate_sigma and not its mean.
Masking via threshold is also possible, but do the magnitudes of the eigenvalues have physical meaning?
I'm new to image processing so any help would be greatly appreciated. I have attached the data for my image also (https://drive.google.com/file/d/1Wdvqtv4rUmLtSNiq-p0M5BwSKxUoW76T/view?usp=sharing).
Thanks!


texts: Does mixing brown noise and white noise create pink noise?

texts:
What happens if you take a white noise generator, run it through a one-pole LPF set to a low cutoff frequency, and then mix the two together? Do you get pink noise?

Say you have this code:

float noise = noiseGenerator();
float pinkNoise = pinkingFilter(noise); //3dB/oct filter
float brownNoise = onePoleLPF(noise); //6dB/oct filter

Would then pinkNoise roughly equal a mix of noise and brownNoise?

Or going further, can you get intermediate colors by just mixing the three noise signals in varying amounts? Like could a mix between pinkNoise and brownNoise create a 4.5dB/oct slope noise?


texts: Characteristic and moment generating function of a random variable interpretation

texts:
I have been studying about moments and cumulants of a random variable. Even though the definitions of characteristic and moments generating function are very similar (only the sign in the exponential is changed) to Fourier and Laplace transforms, respectively.

$\Phi_x(\xi) = E [ e^{j \xi x(\zeta)} ] = \int_{-\infty}^{\infty} f_x(x) e^{j \xi x}  dx $

$\Phi_x(s) = E [ e^{sx(\zeta)} ] = \int_{-\infty}^{\infty} f_x(x) e^{sx}  dx $

The book I've been reading says the variable ξ is not and should not be interpreted as frequency!

How is it supposed to be interpreted? The functions $\Phi_x(\xi)$ and $\Phi_x(s)$ possess all the properties associated with the Fourier
and Laplace transforms, respectively.

Were these exponentials placed there just to be able to perform the Taylor series decomposition and then take the mth-order derivative of every term and extract the appropriate mth-order moment ?


texts: How to extract the co-ordinates of an ECG signal from a scanned image

texts:
I checked the earlier question on "How do I recover the signal from an ECG image"

Currently I am trying to write a program to extract the x and y co-ordinates of an ECG signal embedded within the background grid . I was able to use the image from the earlier post and used some inbuilt functions to crop the image and convert the image to grayscale and then to binary.

This is how I did it:

% Cropping the image and saving the co-ordinates
img=imcrop(imread('ecg_test.jpg'));
% select the boundries in the pop up window , hit on crop image to save the
% cropped image x,y co-ordinates
img1=rgb2gray(img);  % Converting the cropped image to grayscale
Bw=im2bw(img1,0.8);  % converting the grayscale image to binary image
imshow(Bw) ![cropped and binary image of the actual ecg_test.jpg][4]

Now I want to clean the image and then extract the co-ordinates of this image.  How do I implement the MATLAB part of:

img = ColorConvert[
   ImagePad[
    Import["http://i.stack.imgur.com/500Kg.jpg"], {{0, -20}, {0, 0}}],
    "Grayscale"];

Image[
 Transpose[Function[With[{m = Min[#]},
     Map[Function[{v}, If[v == m, 1, 0]], #]]] /@
   Transpose[ImageData[img, "Real"]]
  ]
 ]


texts: Efficient implementation of 2-d circularly symmetric low-pass filter

texts:
How can a two-dimensional ideal circularly symmetric low-pass filter or its approximation be efficiently implemented on data sampled on a square grid? I'm referring to an ideal filter with a spatial frequency response that equals $1$ inside radius $\omega_c$ (the cutoff frequency) and zero outside it, and with ideal impulse response:

$$h[x,y] = \frac{\omega_c}{2\pi \sqrt{x^2 + y^2} } J_1 \big( \omega_c \sqrt{x^2 + y^2} \big).$$

A "not efficient" approach is to convolve with a 2-d filter kernel that is the product of $h[x, y]$ and a circularly symmetric window function.


texts: How to find the instantanous SNR over Fading channel

texts:
I would like to prove that the instantnous SNR of flat fading channel
$$y=hx+n$$
is
$$\gamma=\frac{|h|^2}{N_0}P$$

where $E\{z\}$ is the expectation of $z$ and

$$
E\{x\}=P
$$
$$
E\{n\}=0
$$
$$
E\{n^2\}=N_0
$$

We know that $h$ and $n$ are random varibles so

The noise power is
\begin{align}
E\{n^2\}&=N_0
\end{align}

The problem for me is in the Received signal power

\begin{align}
E\{|hx|^2\}&=E\{hxh^*x^*\}\\
&=E\{|h|^2\}P
\end{align}

Why they take $E\{|h|^2\}=|h|^2$ and we know that $|h|^2$ is random varible and it has it own expectation.

Also what if a noise $n_1$ is given by

$$ n_1=hn+n$$

what is the noise power of $n_1$ where h is the fading and $n$ is nose with zero mean and power $N_0$.

Thanks.


texts: Does taking the Hilbert transform extract the envelope of audio signals?

texts:
I want to extract the envelope of music signals with a sample rate of 44.1kHz.

I used the MATLAB command:  abs(hilbert(mysignals)), but the resulting signal is similar to the original signal, not the envelope.

Does taking the Hilbert transform extract the envelope of an audio signal? Or is my calculation incorrect?


texts:
Prove that with appropriate choices of the phase $\varphi$, one can recover the quadratures of $m(t).$

Attempt

The context of this question is weak optical signals combined with a strong local oscillator at a beamsplitter. And $m(t)$ refers to the complex
modulating signal.

The signal field is given by:  $f(t)=Re \Big[ m(t) \exp (j2 \pi \nu_c) \Big],$ with $m(t)=m_1 (t) + j m_2 (t).$

This is linearly combined with the stronger coherent local oscillator which has the form:

$$c(t) = Re \{ A \exp \Big[ j(2 \pi \nu_{het} t - \varphi) \Big] \},$$

where $\nu_{het} = \nu_c + \Delta,$ has been shifted from the carrier frequency by $\Delta.$ So the output field would be:

$$g_1 (t) = Re \{ R \ A \exp \Big[ j(2 \pi \nu_{het} t - \varphi) \Big] + T \ m(t) \exp (j 2 \pi \nu_c) \},$$

where $R$ and $T$ are reflection and transmission coefficients for the beamsplitter (I was told that for this problem we can assume for simplicity that $R$ is real and $T$ is imaginary).

I am unsure how to proceed with the proof so any explanation would be appreciated.


texts: Normalize Fourier Transform For Video Game

texts:
I am in the process of developing a video game where I would like to sync a character's mouth with an audio clip of their speech that is playing.  I found this great website that shows the profiles of different vowel sounds in a frequency spectrum: http://hyperphysics.phy-astr.gsu.edu/hbase/music/vowel.html

After a little (actually, a lot of) effort, I am able to show in real time the spectrum of the audio that is playing.

The next step I am planning to do is find the 4 largest peaks and try to match it to one of the vowel profiles from the website above.  However, I have noticed that the relative power of the frequency domain coefficients is always strong lower than higher.  The picture I showed is in silence. During an SSSSSS sound, the spectrum does show peaks at higher frequencies:

If I ignore the first few points in the silence FT and then relatively graph the rest, the spectrum has the same shape.  If I ignore a few more, again the lower frequencies are much stronger than the higher frequencies.

I assume that this is the result of the math which I can barely get through to code the basic FFT.  Should I be normalizing the result somehow to be able to see peaks in the higher frequencies?  What process should I use for this?

Thanks in advance for any help you can provide!

To add info requested in the comments, by "in real time", I mean that every frame (at about 60 frames per seconds), I get a set of 1024 samples (sampled at 44100 Hz) that is centered at the point of audio that is currently playing).  I then compute the FT of these samples and draw the spectrum.


texts: Minimal frame sync in GNURadio for packet TX/RX?

texts:
I'm trying to do a minimal packet TX/RX using GFSK.
Signal at RX seems to be OK, but I can't get frame sync to work. This is my flowgraph and the signal I'm seeing:

As you can see, signal/bits at TX and RX are perfectly identical. However, Correlate Access Code block is not getting a match on access code at all.
I found that CAC block is deprecated, but struggling to find a simple alternative.
Example code (https://github.com/gnuradio/gnuradio/blob/master/gr-digital/examples/packet/packet_rx.grc) is quite complex (with FEC support and everything), and I'm not sure if this is the only way to go.
Is there any simpler way to do a frame sync?
I initially started with a Packet Encoder/Decoder (which are also deprecated), and newer way seems to be more and more complex...


texts: Reason of Fast Fourier and Inverse Fast Fourier

texts:
I am not very good at mathematics
I was doing some image processing so I came to know about FFT and IFFT
I was learning about them but the mathematics involved seems to be too dense
So can someone tell me about what are FFT and IFFT and their general application in short in LayMan's Language
So that again when I go to learn it I get it a little easier.
thanks


texts:
I'm trying to make a robot dance on several songs using Arduino.
I need to understand how to get the rhythm from the songs, something like detecting the beats, calculating the bpm..
I'm not expert but I read that probably I'll have to use some filters, is that right? or I can use a simple mic and opamp? do I have to analayze some frequencies?
I hope someone has some ideas and can explaing me how to do such a thing, thanks!! :)


texts: Haralick features, am I approaching this correctly?

texts:
I'm trying to personalise a custom CBIR by adding more features not only based on colour, but based on textures' features.

What I don't know is if I'm approaching this well with Haralick's features.

I'm doing this:

I've used 10 images, 8 for different "textures" and 2 for tests:

All px.jpg are sand images, and all rx.jpg are squares. The test image is squares and test2 image is sand. (NOTE, I don't want to classify, I want a value for each image).

Well, when I query the image test (squares) on the rest of the images (r1, r2, r3, r4, p1, p2, p3, p4) I got the Haralick's features such as Angular 2nd moment, Contrast, Correlation... but I can't find a relationship with the squares, or a not-relationship with sand. I'm making a euclidean distance in this case:

    double features[15]; // I'm introducing other things as well
    Haralick h;
    h.CalculateCoocurrenceMatrix(features);
    correlation();

    for (int j = 1; j<=11; j++) {
        float result = sqrt(pow(featurestest[j]-features[j],2));
        if (result > f_max[j]) {
            f_max[j] = result; // to normalize later
        }

        jsonFeatures["F"+to_string(j)] = result;
    }

My question is: am I doing this approach well? if not, how should I do it?

Thank you in advance.

EDIT:
If it wasn't clear, I don't want to classify such like I'm doing in the example. I'm using "sand" and "squares" only to make the results more visually when sorting (for example seeing p1 p2 r1 p3 p4 r2 p5 r3 r4 r5 on the final match feature).


texts: Triangle-kernel filtering problems at non-integer ratios

texts:
I'm trying to implement image downscaling using a triangle filter (it's narrow therefore fast, and goes hand in hand with linear interpolation when upscaling), and the triangle function used for filtering has a span proportional to the downscaling ratio.

Alas when the ratio is non-integer, grid artifacts appear (as shown above) on the downscaled image. The following plot shows what happens with the way triangles for each pixel of the input image overlap and add up in those cases:

I don't know what I should do to fix/work around this, but I'm sure this problem must have been encountered by others before (Google was no help).


texts:
i am working on a receiver in the space communication. Under low SNR, big frequency offset, phase offset. I started with phase/frequency correction. I used digital PLL and it doesn't work for it. I have read in the book by Rice about "cycle slip" and I can use slip detector in such case. But I dont understand how such detector can be implemented. Can someone provide some information /explanation about it?
Does anyone know the simplest algorithm for cycle slipping correction?
I have found in the book by Rice ( Appendix C) about the cycle slip -effect. I have never heard about it, but it describes perfectly what I get in my simulation of phase correction. Rice didn't provide information about how to correct this effect but mentioned it.


texts:
I am learning Fourier analysis and without any teacher, just trying to read books on my own. I think I have made some decent progress but they are a couple of points which are still very unclear to me and that I can't find explained in any of the books I have.

One of the sources that I found in this document, which is great because it teaches DFT without really using complex numbers. I shall say that I understand complex numbers and I am aware of Euler's Formula.

So what they say in this PDF/document is that in fact, in the simple case you can create a DFT by just using $N/2$ coefficients. I thought the choice of $N/2$ was related to the Nyquist frequency. If the signal contains $N$ samples then the Nyquist sampling theorem says that the signal can't contain a wave whose frequency is higher than half of the sampling frequency (hence the $N/2$ harmonics in the DFT).

So to me, explained that way everything made a lot of sense and in the simplest case you just needed to do something like this:

\begin{align}
a[k] &= \sum_{x = 0}^{N-1} s[x] \cos\left(2 \pi k {1 \over N } x\right), \quad \text{ for } k = \left\{ 1, 2, ..., \frac N2\right\},\\
b[k] &= \sum_{x = 0}^{N-1} s[x] \sin\left(2 \pi k {1 \over N } x\right), \quad \text{ for } k = \left\{ 1, 2, ..., \frac N2\right\}.
\end{align}

So this seemed simple. Now it says that when $k = 0$ and when $k = N/2$ then we need to divide $a$ and $b$ by $N$ or multiply them by $2/N$ otherwise. I understand why when $k = 0$, because it's the DC offset, but didn't really understand why you had to do the same thing when $k = N/2$ until I read this post

QUESTION 1: It seems to indicate that when you use the exponential form of the DFT then when $k = N/2$ then you have $\exp(\pi)$ which is equal to 1. Then it seems that in that situation the coefficient $a[N/2]$ has a particular meaning by I don't know which one?

Now this is where I am lost. In the "complete" equation for the DFT you don't compute $N/2$ coefficients by $N$ coefficients. That means that has soon as $k > N/2$ then the frequency of the harmonics is greater than the Nyquist frequency. I have illustrated this with the following image:

We have $N=8$ samples, thus the fundamental frequency is $\frac 18$ and we have harmonics: $1\cdot \frac 18, \ 2\cdot \frac 18, \ 3\cdot \frac 18, \ 4\cdot \frac 18$. However as soon as go above that then the harmonics go beyond the Nyquist frequency.
QUESTION 2: why do we test the signal with harmonics whose frequency go beyond the Nyquist frequency?

Finally, and I think this is actually related to question 2, I keep reading about positive and negative frequencies, but I just can make sense of this at all? This is my question 3.
QUESTION 3: could you please briefly explained if that's possible why we do speak and need positive and negative frequencies? Why do we need to care for negative frequencies?


texts:
What is the purpose of Balance / Fader in a digital Audio amplifier?
In a signal flow what is the most apt position to place balance and fader(for example after before a compressor unit or at the starting of the signal flow or at the end of signal flow) and why?


texts: Preserving integral through downsampling

texts:
I've got a high-resolution image segmentation, which I'm downsampling (and non-linearly deforming) to fit on an ROI image. I typically blur the high-res image before downsampling with a Gaussian-fwhm related to the downsampling factor. The result is high-res and low-res images with the same mean.

I'd like now for my high-res and low-res images to have the same sum. If I scale up the low-res image by the downsampling factor, will local sums be preserved across the image? Is there a better way to downsample and preserving local sums?


texts: How to reconstruct RF signal using IQ data

texts:
I have a complex IQ signal: $s =  X  + Yi$. I know my original $f_0$ was $5\textrm{ MHz}$.

I can calculate the magnitude/amplitude of $s$: $A = \sqrt{X^2 + Y^2}$
and phase, $\phi = \arctan\left(\frac YX\right) $

But how can I reconstruct the signal using this frequency, amplitude and phase data?

I need the output as a time varying signal which incorporates the amplitude and phase information from the IQ data. I have found this:

  The reconstruction of RF-data from IQ-data is straightforward. It is a
  reversal of the complex demodulation  in  the  previous  section.  The
  decimation  is  reversed  by  interpolation.  The  low- pass  filter
  cannot  be  reversed,  but  should  be  chosen  without  loss  of
  information  in  the  first place.  The  down  mixing  is  reversed
  by  up  mixing.  At  last,  the  RF-signal  is  found  by  taking the
  real-value of the complex up-mixed signal.

But I could really do with a worked example.

EDIT:

This MATLAB code appears to do what I want:

IQ = resample(IQ,downsampl,1);

[Z,X] = size(IQ);
time = [0:Z-1]'*1/F * ones(1,X);

signal = IQ .* exp(+2*pi*i*time*carrier);

signal = real(signal);

where resample resamples IQ at downsampl/1 times.

It is not clear to me what I should use for downsampl and what F is?


texts: About using Walsh matrix as spreading sequence

texts:
I'm asking about some details about using the Walsh matrix as spreading sequence code. For example, suppose I'm using a $4\times 4$ Walsh matrix given by

$$H = \begin{pmatrix}
1 & 1 & 1 & 1 \\
1 & -1 & 1 & -1 \\
1 & 1 & -1 & -1 \\
1 & -1 & -1 & 1
\end{pmatrix}$$
and I have data $x$ which will be spread on that Walsh matrix.

My question: what is the length of $x$ supposed to be? Must it be less than 4? Or can it be more than 4 too? In other words, suppose I have a signal $x = 10$ , it's supposed that is easy to spread the signal $x$ on the first column of $H$. The question is,  What's about if the signal $x = 101101$ , can we do the spread on the first column too? where the signal $x$ has length larger than the first column.


texts:
I am bit unsure whether i am doing i sine sweep correctly.
I am feeding my system a sine wave with f hertz within a given time interval and log its output.

And then i do the same for another frequency.

I tried this technique but seem to do it incorrect or so i think.
I tested for 2 frequencies, and fed them into  Ident in matlab and got a fit rate of 98 % which makes no sense that it could build a such a good bode plot based on 2 values..

I am doing something wrong.. should it not have more than one 2 input or is it enough?


texts: How to transform a Fractional Order Laplace Transfer Function into a digital filter?

texts:
I'm working with loudspeaker impedance analysis. Electrical behavior of loudspeakers can be modeled with RLC networks. But real loudspeakers have components, that exhibit some non-linear and frequency dependent behaviors, that make them difficult to model with simple LTI systems.

One of the problems with loudspeakers is, the voice coil inductance decreases with frequency. (ignoring low frequency behavior, for this question, and focusing on higher freq's). Understanding impedance Z(s) as a Laplace Transfer Function between current and voltage signals, V(s) = Z(s) * I(s).
Speaker impedance does not follow the simple formula of a RL series circuit, Z(s) = Re + (L * s). Neither behaves like a plain resistance Z(s) = Re .

Instead it follows a formula where the s is raised to a non-integer exponent. So it would be Z(s) = Re + (L * (s^a)).

Impedance rises with frequency, not proportionally with f, but proportionally with f^a. Where f = frequency, and a is a real fractional number.

in practice this number called a is around 0.7 (depending on spkr). This phenomenon occurs because frequency dependency of magnetic permeability of the iron core.

following picture describes the issue:

See at 5KHz, impedance is 17 Ohms. At twice this frequency, 10KHz. With a simple inductance, a doubled impedance would be expected, 34 Ohms. But that does not occur. Impedance is increased but not doubled. It is around 27 Ohms. Which is a smaller value than expected with a plain inductance.

Now, what I want to do is to transform the Laplace TF Z = Re + L * (s^0.7), to a discrete-time, Z-transform TF, and then to a IIR digital filter. That would allow to see and analyze the current waveform, from a given voltage signal. Voltage signal is a music MP3 file. Because audio amps outputs are voltage controlled.

With integer exponents in the Laplace TF, is very easy to transform with the Match-Z or Tustin methods. But I have no clue how to do it with fractional exponents. Suppose I want to do Match-Z, how I find the roots? Suppose Tustin Method, replacing s with ((z-1)/(z+1)), How can I raise the ((z-1)/(z+1)) term to a 0.7 exponent???

I know this is a bit hard. Thanks in advance.


texts:
I am trying to match images from a large dataset which exactly match the given input image.(images with deformations/transformations are treated as different)

I have tried euclidean distance but it isn't very scalable. I have also tried matching gradient magnitude and orientation pixel by pixel, but again it seems to work in approximate manner and is slow.

I believe, pixel by pixel matching wouldn't be apt for matching from such large dataset. Moreover, even if certain approximation is allowed, can the threshold be decided dynamically rather than prefixed ?


texts: Rectangle detection in "real life" images

texts:
I have seen a lot of posts here and on SO and on various websites, using Matlab, OpenCV, name it, on how to detect rectangles in an image. I have done it myself too, using Matlab. My usual steps were:

Grayscale the image (if not done)
Blur the image (if needed)
Apply an edge detector (e.g. Canny, LoG, etc., tried a lot of them)
Apply the various steps of Hough Transform to find lines
Use some logic in order to tell which line combinations are forming rectangles

I have no problem with these. The problem is that (and with most of other posts I found) it works with relatively simple images, generally synthetics one with few textures, or with high contrast (like a black rectangle surrounding text on a white sheet of paper). Well, maybe that's only me, but I can't make it work with more complex images, where the edgemap might be quite noisy, where I can't afford to manually tune all sort of parameters and thresholds to "make it work".

Here is a dummy project to illustrate the kind of situation I want to address. Imagine that I have a lot of books (dozens) layered and overlapping on a table and I have a camera above, at the right distance to view all of the scene and only the scene. There is no distortion of lens, no blur from the lighting and the resolution is pretty good. The books may be rotated at any angle, they may be of different sizes (though one could assume there is a min and max size), their cover may or may not be different from each others, with possibly images, text, etc. but no colour information is present.

I want to be able to identify the location of every book in that scene. Then, knowing their location, I will be able to extract the appearance of their covers and maybe retrieve their names from a database using these covers (just an example, this part is not important now). Oh, I cannot turn the problem around and use the database to search for the presence of a given book (it would be quite long and I try to constrain this example to rectangle detection).

Now, imagine that, given all the information known about the problem at hand, I conclude that the best solution would be to search for the presence of rectangles in the scene, knowing that books are indeed rectangles and finding all the rectangle could give me a good start. Oh, and I want to do it in real time.

What would you do? I have looked at a some papers from the literature, a lot are doing contour detection to detect horses, mugs, bottles, but they are quite slow to detect the shapes and more I have noticed that often just plain simple rectangles are, well, too simple, not enough distinctive for the algorithms to find.

I am not looking for a ready made answer. I just wonder if anyone has a thought on this. I have not included a list of papers or implementation I tested partly because this question is quite long enough and also I don't want to orient the answers. Remember, this example is made up and what I really want is to detect rectangles. And, disclosure, eventually in my career I assume I will have to detect rectangles (or other shapes) so I was quite ashamed that I wasn't able to do it on real images with the knowledge I have right now   :S

Thanks


texts:
The Measurement systems and instruments are frequently calibrated.
Do anyone know the reason behind different calibrations.
I have read that the measurement uncertainty increases with every step in the calibration chain..Is it right and why is that?


texts: Why does an fftshift in the time domain give a different result compared to the frequency domain?

texts:
I am trying to speed up some of my code and for physical reasons I need to do an FFT shift to get the FFT data in the standard DFT fourier representation. In Matlab this is performed in the frequency domain, with a data swap.

I had believed that I found a numerical trick, by performing the shift in the time domain, before transformation I was able to save a reasonable amount of computation time.

The problem is that for some reason the results are slightly different. I am curious as to why.

The Matlab code is as follows for a 1K signal

%1khz
f=1000;
t=0:.0001:5;
x=sin(2*pi*f*t);
%Wave
subplot(4,1,1);
plot(x)
%FFT
subplot(4,1,2)
ffted=fft(x);
plot(real(ffted));
%FFT Shifted
subplot(4,1,3)
shifted=fftshift(ffted);
plot(real(shifted));
%FFT Shifted by preshifting
subplot(4,1,4)
s=zeros(1,length(x));%In languages like C++ this is computationally inexpesnive
for i=1:length(x)
    if (bitget(i,1))%odd evens, works for real signal
        s(i)=-1;
    else
        s(i)=1;
    end
end
chirpshifted=fft(s.*x);
plot(real(chirpshifted));
%Difference
A=sum(real(chirpshifted));
B=sum(real(shifted));
A
B


texts: Interference Diversity for management of inter-cell interference in cellular communications

texts:
What is interference diversity? and how is it different in CDMA and OFDM? I have a problem understanding the concept
Reference book: Fundamentals of Wireless Communications, chapter3/4 Author: David Tse


texts: Discrete Harmonics - Why multiplying digital frequency by k does not get next harmonic

texts:
For continuous time $ e^{jk\Omega_0t} $ gives a complete set of orthogonal harmonics for fourier decomposition but for discrete  $ e^{jk\omega_0n} $ does not form a complete set orthogonal basis set to decompose a signal because... ?

No Answer yet on 21 Nov 2018,

If you know please help :(


texts: Subtract a digital filter from another filter

texts: How to subtract a digital filter from another one, if their lengths are different. How to make the length equal.


texts: Peak detection for signals that contain high frequencies

texts:
I have an question about peak detection. My purpose is creating a plot with an impulse response and the decay line at / above it.
What I see till now is that there are two options for finding peaks:

the simplest version is to just compare point with the points beside it
$ f(n-1) , f(n) , f(n+1) $
The other assumes that there are at least 3 or 4 points per half wave (DSP related peak detection)

But my simple self recorded sweeps and transfer function generate just a impulse response with high frequencies. So the peak with option 1 is a peak of the highest sample value but isn't necessarily the highest audio peak.
Also the situation with 3 or 4 samples don't exist.

I'm not looking to a complete answer or work around(however examples are fine). Some links or names from papers would be nice.


texts:
I'm trying to implement an equal power crossfade between two uncorrelated signals ($A$ and $B$). From what I understand I need to setup the crossfade in the following way

Come up with a function $f(x)$ that determines the gain applied to a signal where $x$ is from $0$ to $1$. And $0$ represents being all the way faded to signal $A$,and $1$ represents all the way faded to signal $B$.

$f(x)$ should satisfy:
\begin{cases}\begin{align}
f(x)^2 + f(1-x)^2 &= 1\\
f(1) &= 1\\
f(0) &= 0
\end{align}\end{cases}

So I could use the square root cross-fade $f(x) = \sqrt{x}$:

Or I could use the $\cos$ cross-fade $f(x) = \cos\left(\frac{π}{2}x\right)$:

My question is where does this $f(x)^2 + f(1-x)^2 = 1$ restriction come from? I know that we interpret sound pressure on a logarithmic scale but I don't understand why a quadratic relationship is required between two signals to maintain equal perceived volume.


texts: 8-PSK baseband signal with real and imaginary part

texts:
I want to draw the complex-valued baseband signal with real $S_{BB,I}(t)$ and imaginary part $S_{BB,Q}(t)$ for the following bit sequence: 011111001000

I have looked around and found on wikipedia such a timing diagram but for QPSK:

I understand that the bit sequence of a QPSK is grouped in pairs of two, the first bit of the pair is assigned to I, the second part of the pair is assigned to Q. If I=1, the sine starts at the top. If Q=1, the cosine rises.

How would I need to draw this diagram for an 8-PSK? The bits are grouped in pairs of three. So my bit sequence is 011,111,001,000. I am puzzled how to assign the pairs of three to only two signals I and Q.


texts: issue with the MFCC and GMM for audio recognition

texts:
i am working on a project related to audio events recognition in real time like a door bell, baby crying, footstep, I have 9 categories of sound events
so the first step i did was getting many wav files for the categories i want to recognize (the length of each audio file is 2-3 seconds),then i ran the MFCC function for every single audio file and took the average of it and added it to an excel sheet.

the second step was making a test for that model using the Gaussian methods(by matlab)

here i am making a test by 90%training files and 10% testing files

here i am making a test by 80%training files and 20% testing files

so, why am I getting low recognition results ,and is there anyway to improve it?
thanks


texts: What does "noise power is assumed to be the same at all nodes

texts:
I was reading a paper about communications between nodes, and in the simulation section they gave the condition of the simulation as follow :

"...The noise power is assumed to be the same at all nodes..."

In the simulation it is never stated that the transmit power of all nodes are not the same, but if you put them so that they are equal you get an error. My guess is that they mean the noise variances $σ^2$ of the nodes are the same. However the simulation result is way off the mark from the result in the paper.

Did I understand their meaning properly?


texts: How can we generate an orthogonal sinusodial signal to any other signal?

texts:
I have three signals which are
\begin{align}
A&=\cos\left(2\pi ft + \frac \pi2\right)\\
B&=\cos\left(2\pi ft\right)\\
C&=\cos\left(2\pi ft + \frac {5\pi}4\right)
\end{align}

What i need to do is to draw two orthogonal sinusoidal signal D and E which are also orthogonal to the previous three. That is, D and E will be orthogonal to each other and they will be orthogonal with A, B and C too, at the same time. There is no constraint on the carrier frequency. How can i do it?


texts:
Below response shows sinc1,sinc2 and sinc3 filter response.
I understand that notches in the frequency response happens at the output data rate of the sigma-delta ADC.

But I am not able to figure out why there is a notch ? All articles I read mentions that notch happens at integral multiples of output data rate but doesn't give explanation.
Also why notch width increases as order of the sinc filter increases ?

Please let me know if I am missing basic concepts that I should read before finding answer to this question.


texts: Dominant eigenvectors of an unknown matrix

texts:
Do you have any idea about how we can find the principal eigenvectors of an unknown matrix ${H}$?

The elements of $H$ are unknown in general. If you are familiar with channel estimation procedure in wireless communications, $H$ is channel matrix, and every element $(i,j)$ represents the channel gain between antenna element $i$ of transmitter and antenna element $j$ of the receiver. The following are the information that we have:

we do not know the channel matrix,
we know second order statistics, covariance matrix of the channel,
the channel matrix size tends to infinity (large number of antennas in Tx and Rx),
the channel is sparse, that is, it has only two or three dominant eigenvectors regardless of the channel matrix size.

Note that the power method is not a good solution for me, since its complexity is huge, especially because it requires many iterations to find the dominant eigenvectors.

Some information about the interplay between channel and its covariance:

We consider Rayleigh correlated channel coefficients $h_k\sim CN (0,\mathbf{R}_k)$, mutually independent across the users $k$, and $CN$ is circularly symmetric complex Gaussian. We let $\mathbf{R}_k = \mathbf{U}_k \mathbf{\Lambda}_k \mathbf{U}_{k}^{H}$, where $\mathbf{U}_{k}$ is a tall matrix of the eigenvectors corresponding to the non-zero eigenvalues of $\mathbf{R}_k$, given as the diagonal elements of the diagonal matrix $\mathbf{\Lambda}_{k}$.

The Karhunen-Loeve representation of the user channel vector is given as $\mathbf{h}_{k} = \mathbf{U}_{k}\mathbf{\Lambda}_{k}^{\frac{1}{2}}\mathbf{w}_{k}$, where the entries of $\mathbf{w}_{k}$ are i.i.d. with zero mean Gaussian distributed with variance 1.


texts: Help in interpreting the auto correlation graph

texts:
I want to check if a time series is (a) random  (b) independent. For these I am using the autocorrelation (AC). Autocorrelation
refers to the correlation of a time series with its own past and future values.

If random variables are independent then there should be no correlation between them. The reverse is not true.

I have two different vectors as observations obtained from sensor measurement that are a sequence of samples and I treat this as a time series -denoted by variables phi and phit respectively.

Question 1: I cannot interpret the graph, what the spike indicates. Please help
Question 2: For time series generated from linear system in order to check for randomness, what  is the test to perform and what is the value to expect??
Question 3: For time series obtained from linear system in order to check for independence between samples, what test to perform and what is the value to expect?
Question 4: The given time series is obtained from nonlinear dynamical system that is chaotic.  If correlation test is not valid, then what tests should I perform to check for (a) randomness (b) independence?


texts: Whats is the difference between FIR/IIR filters and Chebyshev/Butterworth filters

texts:
I am new to Signal Processing. From my understanding -- FIR/IIR just refer to the placement of poles and zeros in the z-domain helping us achieve convolution, if FIR and ??? in IIR.

Chebyshev and Butterworth define certain constraints under which a filter can be designed and is indicative of the ripple effect in the pass-band and stop band.

My question is - how are these two related. Can we design an IIR using the Chebyshev constraint. In real-time, is this constraint used using the Impulse Invariance Method or Bilinear Transform for IIR design?

Any answer to this would really help. Thanks


texts:
I'm looking to compute an IDFT but my input is very sparse:
Example:

IDFT length: 7,408,800 (complex floats)
Sparsity: 96.61% to 99.99%

I found this Sparse FFT website but it looks like the library deals with sparse inputs to standard DFT - not IDFT - also it doesn't look like it is not for consumer use.
Are there any techniques I can employ to improve performance for sparse IDFT - before even changing my DFT library? And now, are there any available libraries that give performance improvements for spare IDFT?